mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-24 13:23:41 +01:00
Squashed commit of the following:
commit 96715abd7bc0645b994fc4fff1c764e7bec3042e
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu Sep 6 12:08:34 2018 +0200
Tasks
commit e43b1e2faff5ad1ec4ecd2eb99e8d86e0000a9a2
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu Sep 6 11:47:02 2018 +0200
Fixed print
commit c498c8944b5695b953909711a2e9a542e1bc6015
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu Sep 6 11:11:51 2018 +0200
PT2 and shiftedBk fixed
commit 965cf0361d54df096c9bfca93f9d50ecd946c198
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed Sep 5 18:48:59 2018 +0200
Shifted Bk multistate broken
commit 87ef641b65a122fa1256605cbe098c1a0de04bc0
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed Sep 5 17:23:38 2018 +0200
PT2 fixed
commit a2adb533bcaf96191a24ff5fcef86cd14ac00697
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed Sep 5 16:55:05 2018 +0200
Working on PT2 (broken)
commit 33f52991b65ac42ced5521ef0e713df765268860
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed Sep 5 12:13:23 2018 +0200
Fixed missing argument
commit 712bf75f76421880299dc65b30acfda3d531f709
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed Sep 5 11:42:31 2018 +0200
Fixed floating invalid in PT2
commit cf2412ebd99f9acf573a5180603611f1c3e35155
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed Sep 5 11:34:37 2018 +0200
n_states_diag >= n_states
commit bb415435e4d9be72a285ec789c3d63046a9173e4
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed Sep 5 11:23:41 2018 +0200
Fixed final print
commit f2d339cf7b
Merge: cfa8e1dc 52ca18c1
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 20:51:14 2018 +0200
Merge branch 'thesis' of ssh://github.com/garniron/quantum_package into thesis
commit cfa8e1dc34
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 20:50:09 2018 +0200
restored relative_error
commit bccad69c77
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 20:47:11 2018 +0200
uninitialized variable
commit 52ca18c152
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 20:47:11 2018 +0200
uninitialized variable
commit ba0094f5f8
Merge: 093e3fd0 68458296
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 20:07:15 2018 +0200
merge with garniroy
commit 68458296dc
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue Sep 4 18:43:39 2018 +0200
Almost working but still broken
commit 9ebb88cbf3
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue Sep 4 18:05:00 2018 +0200
Cleaning
commit 873035e016
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue Sep 4 17:31:45 2018 +0200
Squashed commit of the following:
commit 4b9c435dce0f3b3078d573e66fd32b40fca26497
Merge: 74e559c8 093e3fd0
Author: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue Sep 4 16:58:51 2018 +0200
Merge branch 'thesis' of git://github.com/garniron/quantum_package into garniron-thesis
commit 093e3fd021
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 16:13:00 2018 +0200
removed ungodly hack
commit 8529a0f3f6
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:57:19 2018 +0200
reduced prints in pt2_stoch
commit 03b8f353bd
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:41:46 2018 +0200
teeth building check for pt2_stoch
commit 0d91b9310a
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:35:04 2018 +0200
timestamp of first pull
commit 34d9fa0165
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:27:10 2018 +0200
potential numerical precision bug
commit 9a0f900d8c
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:09:51 2018 +0200
tests if teeth can be built
commit dda0dc34df
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 17:48:04 2018 +0200
corrected pt2_find_sample
commit a521f0cb82
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 16:08:02 2018 +0200
tasks get by batches of Nproc
commit 997a5a1265
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 14:18:04 2018 +0200
buffered task_id send
commit 99ea7948e0
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 12:29:12 2018 +0200
unbalanced fragmentation
commit abb3b7e08b
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sun Sep 2 17:18:44 2018 +0200
overflow of pt2_J
commit 8df49f394b
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sun Sep 2 15:58:48 2018 +0200
removed useless computation of intermediate checkpoints
commit 4ba5b79eb3
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sun Sep 2 15:50:14 2018 +0200
dressing only sent for chosen checkpoint
commit a4a6a69459
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sat Sep 1 17:01:56 2018 +0200
cumulative dot_F
commit 6a7f04cb79
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sat Sep 1 16:58:07 2018 +0200
simpler purge
commit 168ca2f2e2
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 21:07:01 2018 +0200
task list optimized
commit de4a0d0caf
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 18:57:03 2018 +0200
removed print
commit fee31d4e3e
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 18:56:23 2018 +0200
dress fragmentation
commit 02893a419d
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 15:52:16 2018 +0200
bug in blocked search - replaced with thesis version
commit bb6e073cf1
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 21:24:45 2018 +0200
ungodly hack to prevent double providing
commit 0609e8c627
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 20:52:05 2018 +0200
debugging
commit a254fdd7cf
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 15:24:07 2018 +0200
parallel bug
commit 2a6c1941d4
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 11:43:11 2018 +0200
corrected when relative_error=0d0
commit bac039bdf1
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 10:58:17 2018 +0200
relative error 1d-5
commit aae9d203ec
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 10:07:02 2018 +0200
potential fragmentation bug
commit ad69f39f99
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Wed Aug 29 20:54:58 2018 +0200
dress_zmq re-implemented
commit d78f64732a
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Wed Aug 29 11:30:19 2018 +0200
pt2_stoch re-implemented
commit 4b9b54e19a
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Aug 28 10:24:38 2018 +0200
removed test for phase_mask_bit
commit 3abccca5e3
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 3 23:44:05 2018 +0200
phasemask_bit
commit 093e3fd021
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 16:13:00 2018 +0200
removed ungodly hack
commit 8529a0f3f6
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:57:19 2018 +0200
reduced prints in pt2_stoch
commit 03b8f353bd
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:41:46 2018 +0200
teeth building check for pt2_stoch
commit 0d91b9310a
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:35:04 2018 +0200
timestamp of first pull
commit 34d9fa0165
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:27:10 2018 +0200
potential numerical precision bug
commit 9a0f900d8c
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Sep 4 14:09:51 2018 +0200
tests if teeth can be built
commit dda0dc34df
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 17:48:04 2018 +0200
corrected pt2_find_sample
commit a521f0cb82
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 16:08:02 2018 +0200
tasks get by batches of Nproc
commit 997a5a1265
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 14:18:04 2018 +0200
buffered task_id send
commit 99ea7948e0
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Mon Sep 3 12:29:12 2018 +0200
unbalanced fragmentation
commit abb3b7e08b
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sun Sep 2 17:18:44 2018 +0200
overflow of pt2_J
commit 8df49f394b
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sun Sep 2 15:58:48 2018 +0200
removed useless computation of intermediate checkpoints
commit 4ba5b79eb3
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sun Sep 2 15:50:14 2018 +0200
dressing only sent for chosen checkpoint
commit a4a6a69459
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sat Sep 1 17:01:56 2018 +0200
cumulative dot_F
commit 6a7f04cb79
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Sat Sep 1 16:58:07 2018 +0200
simpler purge
commit 168ca2f2e2
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 21:07:01 2018 +0200
task list optimized
commit de4a0d0caf
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 18:57:03 2018 +0200
removed print
commit fee31d4e3e
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 18:56:23 2018 +0200
dress fragmentation
commit 02893a419d
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 31 15:52:16 2018 +0200
bug in blocked search - replaced with thesis version
commit bb6e073cf1
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 21:24:45 2018 +0200
ungodly hack to prevent double providing
commit 0609e8c627
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 20:52:05 2018 +0200
debugging
commit a254fdd7cf
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 15:24:07 2018 +0200
parallel bug
commit 2a6c1941d4
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 11:43:11 2018 +0200
corrected when relative_error=0d0
commit bac039bdf1
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 10:58:17 2018 +0200
relative error 1d-5
commit aae9d203ec
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Thu Aug 30 10:07:02 2018 +0200
potential fragmentation bug
commit ad69f39f99
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Wed Aug 29 20:54:58 2018 +0200
dress_zmq re-implemented
commit d78f64732a
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Wed Aug 29 11:30:19 2018 +0200
pt2_stoch re-implemented
commit 4b9b54e19a
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Tue Aug 28 10:24:38 2018 +0200
removed test for phase_mask_bit
commit 3abccca5e3
Author: Yann Garniron <yann.garniron@yahoo.fr>
Date: Fri Aug 3 23:44:05 2018 +0200
phasemask_bit
This commit is contained in:
parent
74e559c85a
commit
4098b05202
@ -10,12 +10,11 @@ program fci_zmq
|
||||
|
||||
double precision :: hf_energy_ref
|
||||
logical :: has
|
||||
double precision :: relative_error, absolute_error
|
||||
double precision :: relative_error
|
||||
integer :: N_states_p
|
||||
character*(512) :: fmt
|
||||
|
||||
relative_error=PT2_relative_error
|
||||
absolute_error=PT2_absolute_error
|
||||
|
||||
pt2 = -huge(1.e0)
|
||||
threshold_davidson_in = threshold_davidson
|
||||
@ -72,7 +71,7 @@ program fci_zmq
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1.d0
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
call ZMQ_pt2(CI_energy, pt2,relative_error,absolute_error,error) ! Stochastic PT2
|
||||
call ZMQ_pt2(CI_energy, pt2,relative_error,error) ! Stochastic PT2
|
||||
threshold_selectors = threshold_selectors_save
|
||||
threshold_generators = threshold_generators_save
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
@ -184,7 +183,7 @@ program fci_zmq
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1d0
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
call ZMQ_pt2(CI_energy, pt2,relative_error,absolute_error,error) ! Stochastic PT2
|
||||
call ZMQ_pt2(CI_energy, pt2,relative_error,error) ! Stochastic PT2
|
||||
threshold_selectors = threshold_selectors_save
|
||||
threshold_generators = threshold_generators_save
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
|
||||
@ -1,246 +0,0 @@
|
||||
program fci_zmq
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
double precision, allocatable :: pt2(:)
|
||||
integer :: degree
|
||||
integer :: n_det_before, to_select
|
||||
double precision :: threshold_davidson_in
|
||||
|
||||
allocate (pt2(N_states))
|
||||
|
||||
double precision :: hf_energy_ref
|
||||
logical :: has
|
||||
double precision :: relative_error, absolute_error
|
||||
integer :: N_states_p
|
||||
character*(512) :: fmt
|
||||
|
||||
relative_error=PT2_relative_error
|
||||
absolute_error=PT2_absolute_error
|
||||
|
||||
pt2 = -huge(1.e0)
|
||||
threshold_davidson_in = threshold_davidson
|
||||
threshold_davidson = threshold_davidson_in * 100.d0
|
||||
SOFT_TOUCH threshold_davidson
|
||||
|
||||
call diagonalize_CI
|
||||
call save_wavefunction
|
||||
|
||||
call ezfio_has_hartree_fock_energy(has)
|
||||
if (has) then
|
||||
call ezfio_get_hartree_fock_energy(hf_energy_ref)
|
||||
else
|
||||
hf_energy_ref = ref_bitmask_energy
|
||||
endif
|
||||
|
||||
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
|
||||
N_states_p = min(N_det,N_states)
|
||||
endif
|
||||
|
||||
n_det_before = 0
|
||||
|
||||
character*(8) :: pt2_string
|
||||
double precision :: correlation_energy_ratio
|
||||
double precision :: threshold_selectors_save, threshold_generators_save
|
||||
threshold_selectors_save = threshold_selectors
|
||||
threshold_generators_save = threshold_generators
|
||||
double precision :: error(N_states)
|
||||
|
||||
correlation_energy_ratio = 0.d0
|
||||
|
||||
if (.True.) then ! Avoid pre-calculation of CI_energy
|
||||
do while ( &
|
||||
(N_det < N_det_max) .and. &
|
||||
(maxval(abs(pt2(1:N_states))) > pt2_max) .and. &
|
||||
(correlation_energy_ratio <= correlation_energy_ratio_max) &
|
||||
)
|
||||
write(*,'(A)') '--------------------------------------------------------------------------------'
|
||||
|
||||
|
||||
if (do_pt2) then
|
||||
pt2_string = ' '
|
||||
pt2 = 0.d0
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1d0
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
call ZMQ_pt2(CI_energy, pt2,relative_error,absolute_error,error) ! Stochastic PT2
|
||||
threshold_selectors = threshold_selectors_save
|
||||
threshold_generators = threshold_generators_save
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
else
|
||||
pt2_string = '(approx)'
|
||||
endif
|
||||
|
||||
|
||||
correlation_energy_ratio = (CI_energy(1) - hf_energy_ref) / &
|
||||
(CI_energy(1) + pt2(1) - hf_energy_ref)
|
||||
correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
|
||||
|
||||
N_states_p = min(N_det,N_states)
|
||||
|
||||
print *, ''
|
||||
print '(A,I12)', 'Summary at N_det = ', N_det
|
||||
print '(A)', '-----------------------------------'
|
||||
print *, ''
|
||||
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
|
||||
print *, ''
|
||||
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))'
|
||||
write(*,fmt) ('State',k, k=1,N_states_p)
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))'
|
||||
write(*,fmt) '# E ', CI_energy(1:N_states_p)
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', CI_energy(1:N_states_p)-CI_energy(1)
|
||||
write(*,fmt) '# Excit. (eV)', (CI_energy(1:N_states_p)-CI_energy(1))*27.211396641308d0
|
||||
endif
|
||||
write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
|
||||
write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
|
||||
write(*,'(A)') '#'
|
||||
write(*,fmt) '# E+PT2 ', (CI_energy(k)+pt2(k),error(k), k=1,N_states_p)
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', ( (CI_energy(k)+pt2(k)-CI_energy(1)-pt2(1)), &
|
||||
dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
|
||||
write(*,fmt) '# Excit. (eV)', ( (CI_energy(k)+pt2(k)-CI_energy(1)-pt2(1))*27.211396641308d0, &
|
||||
dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
|
||||
endif
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
print *, ''
|
||||
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
print*, 'correlation_ratio = ', correlation_energy_ratio
|
||||
|
||||
do k=1, N_states_p
|
||||
print*,'State ',k
|
||||
print *, 'PT2 = ', pt2(k)
|
||||
print *, 'E = ', CI_energy(k)
|
||||
print *, 'E+PT2'//pt2_string//' = ', CI_energy(k)+pt2(k), ' +/- ', error(k)
|
||||
enddo
|
||||
|
||||
print *, '-----'
|
||||
if(N_states.gt.1)then
|
||||
print *, 'Variational Energy difference (au | eV)'
|
||||
do i=2, N_states_p
|
||||
print*,'Delta E = ', (CI_energy(i) - CI_energy(1)), &
|
||||
(CI_energy(i) - CI_energy(1)) * 27.211396641308d0
|
||||
enddo
|
||||
print *, '-----'
|
||||
print*, 'Variational + perturbative Energy difference (au | eV)'
|
||||
do i=2, N_states_p
|
||||
print*,'Delta E = ', (CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1))), &
|
||||
(CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1))) * 27.211396641308d0
|
||||
enddo
|
||||
endif
|
||||
call ezfio_set_full_ci_zmq_energy_pt2(CI_energy(1)+pt2(1))
|
||||
call dump_fci_iterations_value(N_det,CI_energy,pt2)
|
||||
|
||||
n_det_before = N_det
|
||||
if (s2_eig) then
|
||||
to_select = N_det/2+1
|
||||
to_select = max(N_det/2+1, to_select)
|
||||
to_select = min(to_select, N_det_max-n_det_before)
|
||||
else
|
||||
to_select = N_det
|
||||
to_select = max(N_det, to_select)
|
||||
to_select = min(to_select, N_det_max-n_det_before)
|
||||
endif
|
||||
call save_natural_mos
|
||||
call map_deinit(mo_integrals_map)
|
||||
FREE mo_integrals_map
|
||||
PROVIDE mo_integrals_map
|
||||
call four_index_transform_block(ao_integrals_map,mo_integrals_map, &
|
||||
mo_coef, size(mo_coef,1), &
|
||||
1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num, &
|
||||
1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num)
|
||||
|
||||
call ZMQ_selection(to_select, pt2)
|
||||
|
||||
PROVIDE psi_coef
|
||||
PROVIDE psi_det
|
||||
PROVIDE psi_det_sorted
|
||||
|
||||
if (N_det >= N_det_max) then
|
||||
threshold_davidson = threshold_davidson_in
|
||||
end if
|
||||
call diagonalize_CI
|
||||
call save_wavefunction
|
||||
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
|
||||
enddo
|
||||
endif
|
||||
|
||||
if (N_det < N_det_max) then
|
||||
threshold_davidson = threshold_davidson_in
|
||||
call diagonalize_CI
|
||||
call save_wavefunction
|
||||
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
|
||||
call ezfio_set_full_ci_zmq_energy_pt2(CI_energy(1)+pt2(1))
|
||||
endif
|
||||
|
||||
if (do_pt2) then
|
||||
pt2 = 0.d0
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1d0
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
call ZMQ_pt2(CI_energy, pt2,relative_error,absolute_error,error) ! Stochastic PT2
|
||||
threshold_selectors = threshold_selectors_save
|
||||
threshold_generators = threshold_generators_save
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
call ezfio_set_full_ci_zmq_energy(CI_energy(1))
|
||||
call ezfio_set_full_ci_zmq_energy_pt2(CI_energy(1)+pt2(1))
|
||||
endif
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
print*, 'correlation_ratio = ', correlation_energy_ratio
|
||||
|
||||
|
||||
call dump_fci_iterations_value(N_det,CI_energy,pt2)
|
||||
|
||||
print *, ''
|
||||
print '(A,I12)', 'Summary at N_det = ', N_det
|
||||
print '(A)', '-----------------------------------'
|
||||
print *, ''
|
||||
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
|
||||
print *, ''
|
||||
|
||||
|
||||
N_states_p = min(N_det,N_states)
|
||||
print *, ''
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))'
|
||||
write(*,fmt) ('State',k, k=1,N_states_p)
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))'
|
||||
write(*,fmt) '# E ', CI_energy(1:N_states_p)
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', CI_energy(1:N_states_p)-CI_energy(1)
|
||||
write(*,fmt) '# Excit. (eV)', (CI_energy(1:N_states_p)-CI_energy(1))*27.211396641308d0
|
||||
endif
|
||||
write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
|
||||
write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
|
||||
write(*,'(A)') '#'
|
||||
write(*,fmt) '# E+PT2 ', (CI_energy(k)+pt2(k),error(k), k=1,N_states_p)
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', ( (CI_energy(k)+pt2(k)-CI_energy(1)-pt2(1)), &
|
||||
dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
|
||||
write(*,fmt) '# Excit. (eV)', ( (CI_energy(k)+pt2(k)-CI_energy(1)-pt2(1))*27.211396641308d0, &
|
||||
dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
|
||||
endif
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
print *, ''
|
||||
|
||||
|
||||
|
||||
end
|
||||
@ -12,30 +12,30 @@ subroutine run
|
||||
integer :: i,j,k
|
||||
logical, external :: detEq
|
||||
|
||||
double precision, allocatable :: pt2(:)
|
||||
double precision :: pt2(N_states)
|
||||
integer :: degree
|
||||
integer :: n_det_before, to_select
|
||||
double precision :: threshold_davidson_in
|
||||
|
||||
double precision :: E_CI_before, relative_error, absolute_error, eqt
|
||||
double precision :: E_CI_before(N_states), relative_error, error(N_states)
|
||||
|
||||
allocate (pt2(N_states))
|
||||
pt2(:) = 0.d0
|
||||
|
||||
E_CI_before = psi_energy(1) + nuclear_repulsion
|
||||
E_CI_before(:) = psi_energy(:) + nuclear_repulsion
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1.d0
|
||||
relative_error=PT2_relative_error
|
||||
absolute_error=PT2_absolute_error
|
||||
|
||||
call ZMQ_pt2(E_CI_before, pt2, relative_error, absolute_error, eqt)
|
||||
print *, 'Final step'
|
||||
call ZMQ_pt2(E_CI_before, pt2, relative_error, error)
|
||||
do k=1,N_states
|
||||
print *, 'State ', k
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'PT2 = ', pt2
|
||||
print *, 'E = ', E_CI_before
|
||||
print *, 'E+PT2 = ', E_CI_before+pt2, ' +/- ', eqt
|
||||
print *, 'E = ', E_CI_before(k)
|
||||
print *, 'E+PT2 = ', E_CI_before(k)+pt2(k), ' +/- ', error(k)
|
||||
print *, '-----'
|
||||
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before+pt2(1))
|
||||
enddo
|
||||
call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1))
|
||||
end
|
||||
|
||||
|
||||
|
||||
@ -1,32 +1,108 @@
|
||||
BEGIN_PROVIDER [ integer, fragment_first ]
|
||||
BEGIN_PROVIDER [ integer, pt2_stoch_istate ]
|
||||
implicit none
|
||||
fragment_first = first_det_of_teeth(1)
|
||||
BEGIN_DOC
|
||||
! State for stochatsic PT2
|
||||
END_DOC
|
||||
pt2_stoch_istate = 1
|
||||
END_PROVIDER
|
||||
|
||||
subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
|
||||
|
||||
BEGIN_PROVIDER [ integer, pt2_N_teeth ]
|
||||
&BEGIN_PROVIDER [ integer, pt2_minDetInFirstTeeth ]
|
||||
&BEGIN_PROVIDER [ integer, pt2_n_tasks_max ]
|
||||
&BEGIN_PROVIDER [ integer, pt2_F, (N_det_generators) ]
|
||||
implicit none
|
||||
logical, external :: testTeethBuilding
|
||||
integer :: i
|
||||
pt2_F(:) = 1
|
||||
pt2_n_tasks_max = N_det_generators/100 + 1
|
||||
do i=1,N_det_generators
|
||||
if (maxval(dabs(psi_coef_sorted_gen(i,:))) > 0.005d0) then
|
||||
pt2_F(i) = max(1,min( ((elec_alpha_num-n_core_orb)**2)/4, pt2_n_tasks_max))
|
||||
endif
|
||||
enddo
|
||||
|
||||
if(N_det_generators < 1024) then
|
||||
pt2_minDetInFirstTeeth = 1
|
||||
pt2_N_teeth = 1
|
||||
else
|
||||
pt2_minDetInFirstTeeth = min(5, N_det_generators)
|
||||
do pt2_N_teeth=100,2,-1
|
||||
if(testTeethBuilding(pt2_minDetInFirstTeeth, pt2_N_teeth)) exit
|
||||
end do
|
||||
end if
|
||||
call write_int(6,pt2_N_teeth,'Number of comb teeth')
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
logical function testTeethBuilding(minF, N)
|
||||
implicit none
|
||||
integer, intent(in) :: minF, N
|
||||
integer :: n0, i
|
||||
double precision :: u0, Wt, r
|
||||
|
||||
double precision, allocatable :: tilde_w(:), tilde_cW(:)
|
||||
integer, external :: dress_find_sample
|
||||
|
||||
allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators))
|
||||
|
||||
do i=1,N_det_generators
|
||||
tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate)**2 + 1.d-20
|
||||
enddo
|
||||
|
||||
double precision :: norm
|
||||
norm = 0.d0
|
||||
do i=N_det_generators,1,-1
|
||||
norm += tilde_w(i)
|
||||
enddo
|
||||
|
||||
tilde_w(:) = tilde_w(:) / norm
|
||||
|
||||
tilde_cW(0) = -1.d0
|
||||
do i=1,N_det_generators
|
||||
tilde_cW(i) = tilde_cW(i-1) + tilde_w(i)
|
||||
enddo
|
||||
tilde_cW(:) = tilde_cW(:) + 1.d0
|
||||
|
||||
n0 = 0
|
||||
testTeethBuilding = .false.
|
||||
do
|
||||
u0 = tilde_cW(n0)
|
||||
r = tilde_cW(n0 + minF)
|
||||
Wt = (1d0 - u0) / dble(N)
|
||||
if (dabs(Wt) <= 1.d-3) then
|
||||
return
|
||||
endif
|
||||
if(Wt >= r - u0) then
|
||||
testTeethBuilding = .true.
|
||||
return
|
||||
end if
|
||||
n0 += 1
|
||||
if(N_det_generators - n0 < minF * N) then
|
||||
return
|
||||
end if
|
||||
end do
|
||||
stop "exited testTeethBuilding"
|
||||
end function
|
||||
|
||||
|
||||
|
||||
subroutine ZMQ_pt2(E, pt2,relative_error, error)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
|
||||
implicit none
|
||||
|
||||
character(len=64000) :: task
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull
|
||||
type(selection_buffer) :: b
|
||||
integer, external :: omp_get_thread_num
|
||||
double precision, intent(in) :: relative_error, absolute_error, E(N_states)
|
||||
double precision, intent(in) :: relative_error, E(N_states)
|
||||
double precision, intent(out) :: pt2(N_states),error(N_states)
|
||||
|
||||
|
||||
double precision, allocatable :: pt2_detail(:,:), comb(:)
|
||||
logical, allocatable :: computed(:)
|
||||
integer, allocatable :: tbc(:)
|
||||
integer :: i, j, k, Ncomb, i_generator_end
|
||||
integer, external :: pt2_find
|
||||
integer :: i
|
||||
|
||||
double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
|
||||
double precision, external :: omp_get_wtime
|
||||
double precision :: state_average_weight_save(N_states), w(N_states)
|
||||
double precision :: time
|
||||
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
|
||||
|
||||
if (N_det < max(10,N_states)) then
|
||||
@ -41,25 +117,8 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
|
||||
state_average_weight(pt2_stoch_istate) = 1.d0
|
||||
TOUCH state_average_weight pt2_stoch_istate
|
||||
|
||||
allocate(pt2_detail(N_states,N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
|
||||
sumabove = 0d0
|
||||
sum2above = 0d0
|
||||
Nabove = 0d0
|
||||
provide nproc pt2_F mo_bielec_integrals_in_map mo_mono_elec_integral pt2_w psi_selectors
|
||||
|
||||
provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors
|
||||
|
||||
computed = .false.
|
||||
|
||||
tbc(0) = first_det_of_comb - 1
|
||||
do i=1, tbc(0)
|
||||
tbc(i) = i
|
||||
computed(i) = .true.
|
||||
end do
|
||||
|
||||
Ncomb=size(comb)
|
||||
call get_carlo_workbatch(computed, comb, Ncomb, tbc)
|
||||
|
||||
pt2_detail = 0d0
|
||||
print *, '========== ================= ================= ================='
|
||||
print *, ' Samples Energy Stat. Error Seconds '
|
||||
print *, '========== ================= ================= ================='
|
||||
@ -97,26 +156,15 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
|
||||
endif
|
||||
|
||||
|
||||
call create_selection_buffer(1, 1*2, b)
|
||||
|
||||
integer :: ipos
|
||||
ipos=1
|
||||
|
||||
integer, external :: add_task_to_taskserver
|
||||
|
||||
do i=1,tbc(0)
|
||||
if(tbc(i) > fragment_first) then
|
||||
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') 0, tbc(i)
|
||||
ipos += 20
|
||||
if (ipos > 63980) then
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
|
||||
stop 'Unable to add task to task server'
|
||||
endif
|
||||
character(len=64000) :: task
|
||||
integer :: j,k,ipos
|
||||
ipos=1
|
||||
endif
|
||||
else
|
||||
do j=1,fragment_count
|
||||
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, tbc(i)
|
||||
task = ' '
|
||||
|
||||
do i= 1, N_det_generators
|
||||
do j=1,pt2_F(pt2_J(i))
|
||||
write(task(ipos:ipos+20),'(I9,1X,I9,''|'')') j, pt2_J(i)
|
||||
ipos += 20
|
||||
if (ipos > 63980) then
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
|
||||
@ -125,7 +173,6 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
|
||||
ipos=1
|
||||
endif
|
||||
end do
|
||||
end if
|
||||
enddo
|
||||
if (ipos > 1) then
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
|
||||
@ -149,23 +196,25 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
|
||||
nproc_target = min(nproc_target,nproc)
|
||||
endif
|
||||
|
||||
call omp_set_nested(.true.)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(shared) NUM_THREADS(nproc_target+1) &
|
||||
!$OMP PRIVATE(i)
|
||||
i = omp_get_thread_num()
|
||||
if (i==0) then
|
||||
call pt2_collector(zmq_socket_pull,E(pt2_stoch_istate), b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, absolute_error, w, error)
|
||||
call pt2_collector(zmq_socket_pull, E(pt2_stoch_istate),relative_error, w, error)
|
||||
pt2(pt2_stoch_istate) = w(pt2_stoch_istate)
|
||||
else
|
||||
call pt2_slave_inproc(i)
|
||||
endif
|
||||
!$OMP END PARALLEL
|
||||
call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2')
|
||||
call delete_selection_buffer(b)
|
||||
|
||||
print *, '========== ================= ================= ================='
|
||||
|
||||
deallocate(pt2_detail, comb, computed, tbc)
|
||||
enddo
|
||||
! call omp_set_nested(.false.)
|
||||
|
||||
FREE pt2_stoch_istate
|
||||
state_average_weight(:) = state_average_weight_save(:)
|
||||
TOUCH state_average_weight
|
||||
@ -177,34 +226,6 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error)
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, computed, sumabove, sum2above, Nabove)
|
||||
integer, intent(in) :: tbc(0:size_tbc), Ncomb
|
||||
logical, intent(in) :: computed(N_det_generators)
|
||||
double precision, intent(in) :: comb(Ncomb), pt2_detail(N_states,N_det_generators)
|
||||
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
|
||||
integer :: i, dets(comb_teeth)
|
||||
double precision :: myVal, myVal2
|
||||
|
||||
mainLoop : do i=1,Ncomb
|
||||
call get_comb(comb(i), dets, comb_teeth)
|
||||
do j=1,comb_teeth
|
||||
if(.not.(computed(dets(j)))) then
|
||||
exit mainLoop
|
||||
end if
|
||||
end do
|
||||
|
||||
myVal = 0d0
|
||||
myVal2 = 0d0
|
||||
do j=comb_teeth,1,-1
|
||||
myVal += pt2_detail(pt2_stoch_istate,dets(j)) * pt2_weight_inv(dets(j)) * comb_step
|
||||
sumabove(j) += myVal
|
||||
sum2above(j) += myVal*myVal
|
||||
Nabove(j) += 1
|
||||
end do
|
||||
end do mainLoop
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine pt2_slave_inproc(i)
|
||||
implicit none
|
||||
integer, intent(in) :: i
|
||||
@ -212,411 +233,289 @@ subroutine pt2_slave_inproc(i)
|
||||
call run_pt2_slave(1,i,pt2_e0_denominator)
|
||||
end
|
||||
|
||||
subroutine pt2_collector(zmq_socket_pull, E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, sum2above, Nabove, relative_error, absolute_error, pt2,error)
|
||||
|
||||
subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2, error)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
|
||||
integer, intent(in) :: Ncomb
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
|
||||
double precision, intent(inout) :: pt2_detail(N_states, N_det_generators)
|
||||
double precision, intent(in) :: comb(Ncomb), relative_error, absolute_error, E
|
||||
logical, intent(inout) :: computed(N_det_generators)
|
||||
integer, intent(in) :: tbc(0:size_tbc)
|
||||
double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth)
|
||||
double precision, intent(in) :: relative_error, E
|
||||
double precision, intent(out) :: pt2(N_states), error(N_states)
|
||||
|
||||
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
double precision, allocatable :: pt2_mwen(:,:)
|
||||
double precision, allocatable :: eI(:,:), eI_task(:,:), S(:), S2(:)
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
integer, external :: zmq_delete_tasks
|
||||
integer, external :: zmq_abort
|
||||
integer, external :: pt2_find_sample
|
||||
|
||||
|
||||
integer :: msg_size, rc, more
|
||||
integer :: acc, i, j, robin, N, n_tasks
|
||||
double precision, allocatable :: val(:)
|
||||
integer(bit_kind), allocatable :: det(:,:,:)
|
||||
integer :: more, n, i, p, c, t, n_tasks, U
|
||||
integer, allocatable :: task_id(:)
|
||||
integer, allocatable :: index(:)
|
||||
double precision :: time0
|
||||
double precision :: time, timeLast, Nabove_old
|
||||
|
||||
double precision, external :: omp_get_wtime
|
||||
integer :: tooth, firstTBDcomb, orgTBDcomb, n_tasks_max
|
||||
integer, allocatable :: parts_to_get(:)
|
||||
logical, allocatable :: actually_computed(:)
|
||||
double precision :: eqt
|
||||
character*(512) :: task
|
||||
Nabove_old = -1.d0
|
||||
n_tasks_max = N_det_generators/100+1
|
||||
double precision :: v, x, avg, eqt, E0
|
||||
double precision :: time, time0
|
||||
|
||||
allocate(actually_computed(N_det_generators), parts_to_get(N_det_generators), &
|
||||
pt2_mwen(N_states, n_tasks_max) )
|
||||
|
||||
pt2_mwen(1:N_states, 1:n_tasks_max) = 0.d0
|
||||
do i=1,N_det_generators
|
||||
actually_computed(i) = computed(i)
|
||||
enddo
|
||||
|
||||
parts_to_get(:) = 1
|
||||
if(fragment_first > 0) then
|
||||
do i=1,fragment_first
|
||||
parts_to_get(i) = fragment_count
|
||||
enddo
|
||||
endif
|
||||
|
||||
do i=1,tbc(0)
|
||||
actually_computed(tbc(i)) = .false.
|
||||
end do
|
||||
|
||||
orgTBDcomb = int(Nabove(1))
|
||||
firstTBDcomb = 1
|
||||
integer, allocatable :: f(:)
|
||||
logical, allocatable :: d(:)
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
allocate(val(b%N), det(N_int, 2, b%N), task_id(n_tasks_max), index(n_tasks_max))
|
||||
allocate(task_id(pt2_n_tasks_max), index(pt2_n_tasks_max), f(N_det_generators))
|
||||
allocate(d(N_det_generators+1))
|
||||
allocate(eI(N_states, N_det_generators), eI_task(N_states, pt2_n_tasks_max))
|
||||
allocate(S(pt2_N_teeth+1), S2(pt2_N_teeth+1))
|
||||
|
||||
pt2(:) = -huge(1.)
|
||||
S(:) = 0d0
|
||||
S2(:) = 0d0
|
||||
n = 1
|
||||
t = 0
|
||||
U = 0
|
||||
eI(:,:) = 0d0
|
||||
f(:) = pt2_F(:)
|
||||
d(:) = .false.
|
||||
n_tasks = 0
|
||||
E0 = E
|
||||
more = 1
|
||||
call wall_time(time0)
|
||||
timeLast = time0
|
||||
time0 = omp_get_wtime()
|
||||
|
||||
call get_first_tooth(actually_computed, tooth)
|
||||
Nabove_old = Nabove(tooth)
|
||||
|
||||
logical :: loop
|
||||
loop = .True.
|
||||
pullLoop : do while (loop)
|
||||
|
||||
call pull_pt2_results(zmq_socket_pull, index, pt2_mwen, task_id, n_tasks)
|
||||
do i=1,n_tasks
|
||||
pt2_detail(1:N_states, index(i)) += pt2_mwen(1:N_states,i)
|
||||
parts_to_get(index(i)) -= 1
|
||||
if(parts_to_get(index(i)) < 0) then
|
||||
print *, i, index(i), parts_to_get(index(i))
|
||||
print *, parts_to_get
|
||||
stop "PARTS ??"
|
||||
end if
|
||||
if(parts_to_get(index(i)) == 0) actually_computed(index(i)) = .true.
|
||||
do while (n <= N_det_generators)
|
||||
if(f(pt2_J(n)) == 0) then
|
||||
d(pt2_J(n)) = .true.
|
||||
do while(d(U+1))
|
||||
U += 1
|
||||
end do
|
||||
|
||||
integer, external :: zmq_delete_tasks
|
||||
if (zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,n_tasks,more) == -1) then
|
||||
cycle
|
||||
endif
|
||||
if (more == 0) then
|
||||
loop = .False.
|
||||
endif
|
||||
|
||||
call wall_time(time)
|
||||
|
||||
|
||||
if(time - timeLast > 5d0 .or. (.not.loop)) then
|
||||
timeLast = time
|
||||
do i=1, first_det_of_teeth(1)-1
|
||||
if(.not.(actually_computed(i))) then
|
||||
cycle pullLoop
|
||||
end if
|
||||
! Deterministic part
|
||||
do while(t <= pt2_N_teeth)
|
||||
if(U >= pt2_n_0(t+1)) then
|
||||
t=t+1
|
||||
E0 = 0.d0
|
||||
do i=pt2_n_0(t),1,-1
|
||||
E0 += eI(pt2_stoch_istate, i)
|
||||
end do
|
||||
|
||||
integer, external :: zmq_abort
|
||||
double precision :: E0, avg, prop
|
||||
|
||||
call do_carlo(tbc, Ncomb+1-firstTBDcomb, comb(firstTBDcomb), pt2_detail, actually_computed, sumabove, sum2above, Nabove)
|
||||
firstTBDcomb = int(Nabove(1)) - orgTBDcomb + 1
|
||||
call get_first_tooth(actually_computed, tooth)
|
||||
|
||||
if (firstTBDcomb > Ncomb) then
|
||||
if (zmq_abort(zmq_to_qp_run_socket) == -1) then
|
||||
call sleep(1)
|
||||
if (zmq_abort(zmq_to_qp_run_socket) == -1) then
|
||||
print *, irp_here, ': Error in sending abort signal (1)'
|
||||
endif
|
||||
endif
|
||||
! exit pullLoop
|
||||
endif
|
||||
|
||||
E0 = sum(pt2_detail(pt2_stoch_istate,:first_det_of_teeth(tooth)-1))
|
||||
if (tooth <= comb_teeth) then
|
||||
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
|
||||
prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
|
||||
E0 += pt2_detail(pt2_stoch_istate,first_det_of_teeth(tooth)) * prop
|
||||
avg = E0 + (sumabove(tooth) / Nabove(tooth))
|
||||
eqt = sqrt(1d0 / (Nabove(tooth)-1.d0) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
|
||||
else
|
||||
eqt = 0.d0
|
||||
tooth=comb_teeth
|
||||
exit
|
||||
end if
|
||||
call wall_time(time)
|
||||
if ( ((dabs(eqt/avg) < relative_error) .or. (dabs(eqt) < absolute_error)) .and. Nabove(tooth) >= 10.d0) then
|
||||
! Termination
|
||||
end do
|
||||
|
||||
! Add Stochastic part
|
||||
c = pt2_R(n)
|
||||
if(c > 0) then
|
||||
x = 0d0
|
||||
do p=pt2_N_teeth, 1, -1
|
||||
v = pt2_u_0 + pt2_W_T * (pt2_u(c) + dble(p-1))
|
||||
i = pt2_find_sample(v, pt2_cW)
|
||||
x += eI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i)
|
||||
S(p) += x
|
||||
S2(p) += x**2
|
||||
end do
|
||||
avg = E0 + S(t) / dble(c)
|
||||
pt2(pt2_stoch_istate) = avg
|
||||
! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969)
|
||||
if(c > 2) then
|
||||
eqt = dabs((S2(t) / c) - (S(t)/c)**2) ! dabs for numerical stability
|
||||
eqt = sqrt(eqt / (dble(c) - 1.5d0))
|
||||
error(pt2_stoch_istate) = eqt
|
||||
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
|
||||
if(mod(c,10)==0 .or. n==N_det_generators) then
|
||||
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', c, avg+E, eqt, time-time0, ''
|
||||
if( dabs(error(pt2_stoch_istate) / pt2(pt2_stoch_istate)) < relative_error) then
|
||||
if (zmq_abort(zmq_to_qp_run_socket) == -1) then
|
||||
call sleep(1)
|
||||
if (zmq_abort(zmq_to_qp_run_socket) == -1) then
|
||||
print *, irp_here, ': Error in sending abort signal (2)'
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
endif
|
||||
time = omp_get_wtime()
|
||||
end if
|
||||
n += 1
|
||||
else if(more == 0) then
|
||||
exit
|
||||
else
|
||||
if ( (Nabove(tooth) > 2.d0) .and. (Nabove(tooth) > Nabove_old) ) then
|
||||
print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, ''
|
||||
Nabove_old = Nabove(tooth)
|
||||
call pull_pt2_results(zmq_socket_pull, index, eI_task, task_id, n_tasks)
|
||||
if (zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,n_tasks,more) == -1) then
|
||||
stop 'Unable to delete tasks'
|
||||
endif
|
||||
endif
|
||||
end if
|
||||
end do pullLoop
|
||||
|
||||
if(tooth == comb_teeth+1) then
|
||||
pt2(pt2_stoch_istate) = sum(pt2_detail(pt2_stoch_istate,:))
|
||||
error(pt2_stoch_istate) = 0d0
|
||||
else
|
||||
E0 = sum(pt2_detail(pt2_stoch_istate,:first_det_of_teeth(tooth)-1))
|
||||
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
|
||||
prop = prop * pt2_weight_inv(first_det_of_teeth(tooth))
|
||||
E0 += pt2_detail(pt2_stoch_istate,first_det_of_teeth(tooth)) * prop
|
||||
pt2(pt2_stoch_istate) = E0 + (sumabove(tooth) / Nabove(tooth))
|
||||
error(pt2_stoch_istate) = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
|
||||
end if
|
||||
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call sort_selection_buffer(b)
|
||||
end subroutine
|
||||
|
||||
integer function pt2_find(v, w, sze, imin, imax)
|
||||
implicit none
|
||||
integer, intent(in) :: sze, imin, imax
|
||||
double precision, intent(in) :: v, w(sze)
|
||||
integer :: i,l,h
|
||||
integer, parameter :: block=64
|
||||
|
||||
l = imin
|
||||
h = imax-1
|
||||
|
||||
do while(h-l >= block)
|
||||
i = ishft(h+l,-1)
|
||||
if(w(i+1) > v) then
|
||||
h = i-1
|
||||
else
|
||||
l = i+1
|
||||
do i=1,n_tasks
|
||||
eI(:, index(i)) += eI_task(:, i)
|
||||
f(index(i)) -= 1
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
!DIR$ LOOP COUNT (64)
|
||||
do pt2_find=l,h
|
||||
if(w(pt2_find) >= v) then
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
end subroutine
|
||||
|
||||
|
||||
integer function pt2_find_sample(v, w)
|
||||
implicit none
|
||||
double precision, intent(in) :: v, w(0:N_det_generators)
|
||||
integer :: i,l,r
|
||||
|
||||
l = 0
|
||||
r = N_det_generators
|
||||
|
||||
do while(r-l > 1)
|
||||
i = (r+l) / 2
|
||||
if(w(i) < v) then
|
||||
l = i
|
||||
else
|
||||
r = i
|
||||
end if
|
||||
end do
|
||||
i = r
|
||||
do r=i+1,N_det_generators
|
||||
if (w(r) /= w(i)) then
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
pt2_find_sample = r-1
|
||||
end function
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, comb_teeth ]
|
||||
BEGIN_PROVIDER[ integer, pt2_J, (N_det_generators)]
|
||||
&BEGIN_PROVIDER[ double precision, pt2_u, (N_det_generators)]
|
||||
&BEGIN_PROVIDER[ integer, pt2_R, (N_det_generators)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of teeth in the comb
|
||||
END_DOC
|
||||
comb_teeth = min(1+N_det/10,100)
|
||||
integer :: N_c, N_j, U, t, i
|
||||
double precision :: v
|
||||
logical, allocatable :: d(:)
|
||||
integer, external :: pt2_find_sample
|
||||
|
||||
allocate(d(N_det_generators))
|
||||
|
||||
pt2_R(:) = 0
|
||||
N_c = 0
|
||||
N_j = pt2_n_0(1)
|
||||
d(:) = .false.
|
||||
|
||||
do i=1,N_j
|
||||
d(i) = .true.
|
||||
pt2_J(i) = i
|
||||
end do
|
||||
call random_seed(put=(/3211,64,6566,321,65,321,654,65,321,6321,654,65,321,621,654,65,321,65,654,65,321,65/))
|
||||
call RANDOM_NUMBER(pt2_u)
|
||||
call RANDOM_NUMBER(pt2_u)
|
||||
|
||||
|
||||
|
||||
U = 0
|
||||
|
||||
do while(N_j < N_det_generators)
|
||||
!ADD_COMB
|
||||
N_c += 1
|
||||
do t=0, pt2_N_teeth-1
|
||||
v = pt2_u_0 + pt2_W_T * (dble(t) + pt2_u(N_c))
|
||||
i = pt2_find_sample(v, pt2_cW)
|
||||
if(.not. d(i)) then
|
||||
N_j += 1
|
||||
pt2_J(N_j) = i
|
||||
d(i) = .true.
|
||||
end if
|
||||
end do
|
||||
|
||||
pt2_R(N_j) = N_c
|
||||
|
||||
!FILL_TOOTH
|
||||
do while(U < N_det_generators)
|
||||
U += 1
|
||||
if(.not. d(U)) then
|
||||
N_j += 1
|
||||
pt2_J(N_j) = U
|
||||
d(U) = .true.
|
||||
exit;
|
||||
end if
|
||||
end do
|
||||
enddo
|
||||
if(N_det_generators > 1) then
|
||||
pt2_R(N_det_generators-1) = 0
|
||||
pt2_R(N_det_generators) = N_c
|
||||
end if
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
subroutine get_first_tooth(computed, first_teeth)
|
||||
BEGIN_PROVIDER [ double precision, pt2_w, (N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, pt2_cW, (0:N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, pt2_W_T ]
|
||||
&BEGIN_PROVIDER [ double precision, pt2_u_0 ]
|
||||
&BEGIN_PROVIDER [ integer, pt2_n_0, (pt2_N_teeth+1) ]
|
||||
implicit none
|
||||
logical, intent(in) :: computed(N_det_generators)
|
||||
integer, intent(out) :: first_teeth
|
||||
integer :: i, first_det
|
||||
integer :: i, t
|
||||
double precision, allocatable :: tilde_w(:), tilde_cW(:)
|
||||
double precision :: r, tooth_width
|
||||
integer, external :: pt2_find_sample
|
||||
|
||||
first_det = N_det_generators+1+1
|
||||
first_teeth = 1
|
||||
do i=first_det_of_comb, N_det_generators
|
||||
if(.not.(computed(i))) then
|
||||
first_det = i
|
||||
allocate(tilde_w(N_det_generators), tilde_cW(0:N_det_generators))
|
||||
|
||||
tilde_cW(0) = 0d0
|
||||
|
||||
do i=1,N_det_generators
|
||||
tilde_w(i) = psi_coef_sorted_gen(i,pt2_stoch_istate)**2 + 1.d-20
|
||||
enddo
|
||||
|
||||
double precision :: norm
|
||||
norm = 0.d0
|
||||
do i=N_det_generators,1,-1
|
||||
norm += tilde_w(i)
|
||||
enddo
|
||||
|
||||
tilde_w(:) = tilde_w(:) / norm
|
||||
|
||||
tilde_cW(0) = -1.d0
|
||||
do i=1,N_det_generators
|
||||
tilde_cW(i) = tilde_cW(i-1) + tilde_w(i)
|
||||
enddo
|
||||
tilde_cW(:) = tilde_cW(:) + 1.d0
|
||||
|
||||
pt2_n_0(1) = 0
|
||||
do
|
||||
pt2_u_0 = tilde_cW(pt2_n_0(1))
|
||||
r = tilde_cW(pt2_n_0(1) + pt2_minDetInFirstTeeth)
|
||||
pt2_W_T = (1d0 - pt2_u_0) / dble(pt2_N_teeth)
|
||||
if(pt2_W_T >= r - pt2_u_0) then
|
||||
exit
|
||||
end if
|
||||
end do
|
||||
|
||||
do i=comb_teeth+1, 1, -1
|
||||
if(first_det_of_teeth(i) < first_det) then
|
||||
first_teeth = i
|
||||
exit
|
||||
pt2_n_0(1) += 1
|
||||
if(N_det_generators - pt2_n_0(1) < pt2_minDetInFirstTeeth * pt2_N_teeth) then
|
||||
stop "teeth building failed"
|
||||
end if
|
||||
end do
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer*8, size_tbc ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Size of the tbc array
|
||||
END_DOC
|
||||
size_tbc = int((comb_teeth+1),8)*int(N_det_generators,8) + fragment_count*fragment_first
|
||||
END_PROVIDER
|
||||
|
||||
subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
|
||||
implicit none
|
||||
integer, intent(inout) :: Ncomb
|
||||
double precision, intent(out) :: comb(Ncomb)
|
||||
integer, intent(inout) :: tbc(0:size_tbc)
|
||||
logical, intent(inout) :: computed(N_det_generators)
|
||||
integer :: i, j, last_full, dets(comb_teeth)
|
||||
integer :: icount, n
|
||||
integer :: k, l
|
||||
l=first_det_of_comb
|
||||
call RANDOM_NUMBER(comb)
|
||||
do i=1,size(comb)
|
||||
comb(i) = comb(i) * comb_step
|
||||
!DIR$ FORCEINLINE
|
||||
call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
|
||||
Ncomb = i
|
||||
if (tbc(0) == N_det_generators) return
|
||||
do while (computed(l))
|
||||
l=l+1
|
||||
enddo
|
||||
k=tbc(0)+1
|
||||
tbc(k) = l
|
||||
computed(l) = .True.
|
||||
tbc(0) = k
|
||||
do t=2, pt2_N_teeth
|
||||
r = pt2_u_0 + pt2_W_T * dble(t-1)
|
||||
pt2_n_0(t) = pt2_find_sample(r, tilde_cW)
|
||||
end do
|
||||
pt2_n_0(pt2_N_teeth+1) = N_det_generators
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
subroutine get_comb(stato, dets, ct)
|
||||
implicit none
|
||||
integer, intent(in) :: ct
|
||||
double precision, intent(in) :: stato
|
||||
integer, intent(out) :: dets(ct)
|
||||
double precision :: curs
|
||||
integer :: j
|
||||
integer, external :: pt2_find
|
||||
|
||||
curs = 1d0 - stato
|
||||
do j = comb_teeth, 1, -1
|
||||
!DIR$ FORCEINLINE
|
||||
dets(j) = pt2_find(curs, pt2_cweight,size(pt2_cweight), first_det_of_teeth(j), first_det_of_teeth(j+1))
|
||||
curs -= comb_step
|
||||
end do
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine add_comb(comb, computed, tbc, stbc, ct)
|
||||
implicit none
|
||||
integer*8, intent(in) :: stbc
|
||||
integer, intent(in) :: ct
|
||||
double precision, intent(in) :: comb
|
||||
logical, intent(inout) :: computed(N_det_generators)
|
||||
integer, intent(inout) :: tbc(0:stbc)
|
||||
integer :: i, k, l, dets(ct)
|
||||
|
||||
!DIR$ FORCEINLINE
|
||||
call get_comb(comb, dets, ct)
|
||||
|
||||
k=tbc(0)+1
|
||||
do i = 1, ct
|
||||
l = dets(i)
|
||||
if(.not.(computed(l))) then
|
||||
tbc(k) = l
|
||||
k = k+1
|
||||
computed(l) = .true.
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
pt2_w(:pt2_n_0(1)) = tilde_w(:pt2_n_0(1))
|
||||
do t=1, pt2_N_teeth
|
||||
tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t))
|
||||
if (tooth_width == 0.d0) then
|
||||
tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1)))
|
||||
endif
|
||||
ASSERT(tooth_width > 0.d0)
|
||||
do i=pt2_n_0(t)+1, pt2_n_0(t+1)
|
||||
pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width
|
||||
end do
|
||||
end do
|
||||
tbc(0) = k-1
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, pt2_stoch_istate ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! State for stochatsic PT2
|
||||
END_DOC
|
||||
pt2_stoch_istate = 1
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, pt2_weight, (N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, pt2_cweight, (N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, pt2_cweight_cache, (N_det_generators) ]
|
||||
&BEGIN_PROVIDER [ double precision, comb_step ]
|
||||
&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ]
|
||||
&BEGIN_PROVIDER [ integer, first_det_of_comb ]
|
||||
implicit none
|
||||
integer :: i
|
||||
double precision :: norm_left, stato
|
||||
integer, external :: pt2_find
|
||||
|
||||
pt2_weight(1) = psi_coef_generators(1,pt2_stoch_istate)**2
|
||||
pt2_cweight(1) = psi_coef_generators(1,pt2_stoch_istate)**2
|
||||
|
||||
pt2_cW(0) = 0d0
|
||||
do i=1,N_det_generators
|
||||
pt2_weight(i) = psi_coef_generators(i,pt2_stoch_istate)**2
|
||||
pt2_cW(i) = pt2_cW(i-1) + pt2_w(i)
|
||||
end do
|
||||
|
||||
! Important to loop backwards for numerical precision
|
||||
pt2_cweight(N_det_generators) = pt2_weight(N_det_generators)
|
||||
do i=N_det_generators-1,1,-1
|
||||
pt2_cweight(i) = pt2_weight(i) + pt2_cweight(i+1)
|
||||
end do
|
||||
|
||||
do i=1,N_det_generators
|
||||
pt2_weight(i) = pt2_weight(i) / pt2_cweight(1)
|
||||
pt2_cweight(i) = pt2_cweight(i) / pt2_cweight(1)
|
||||
enddo
|
||||
|
||||
do i=1,N_det_generators-1
|
||||
pt2_cweight(i) = 1.d0 - pt2_cweight(i+1)
|
||||
end do
|
||||
pt2_cweight(N_det_generators) = 1.d0
|
||||
|
||||
norm_left = 1d0
|
||||
|
||||
comb_step = 1d0/dfloat(comb_teeth)
|
||||
first_det_of_comb = 1
|
||||
do i=1,N_det_generators
|
||||
if(pt2_weight(i)/norm_left < .5d0*comb_step) then
|
||||
first_det_of_comb = i
|
||||
exit
|
||||
end if
|
||||
norm_left -= pt2_weight(i)
|
||||
end do
|
||||
first_det_of_comb = max(2,first_det_of_comb)
|
||||
call write_int(6, first_det_of_comb-1, 'Size of deterministic set')
|
||||
|
||||
comb_step = (1d0 - pt2_cweight(first_det_of_comb-1)) * comb_step
|
||||
|
||||
stato = 1d0 - comb_step
|
||||
iloc = N_det_generators
|
||||
do i=comb_teeth, 1, -1
|
||||
integer :: iloc
|
||||
iloc = pt2_find(stato, pt2_cweight, N_det_generators, 1, iloc)
|
||||
first_det_of_teeth(i) = iloc
|
||||
stato -= comb_step
|
||||
end do
|
||||
first_det_of_teeth(comb_teeth+1) = N_det_generators + 1
|
||||
first_det_of_teeth(1) = first_det_of_comb
|
||||
if(first_det_of_teeth(1) /= first_det_of_comb) then
|
||||
print *, 'Error in ', irp_here
|
||||
stop "comb provider"
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, pt2_weight_inv, (N_det_generators) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Inverse of pt2_weight array
|
||||
END_DOC
|
||||
integer :: i
|
||||
do i=1,N_det_generators
|
||||
pt2_weight_inv(i) = 1.d0/pt2_weight(i)
|
||||
enddo
|
||||
|
||||
pt2_n_0(pt2_N_teeth+1) = N_det_generators
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -22,12 +22,11 @@ subroutine run_pt2_slave(thread,iproc,energy)
|
||||
logical :: done
|
||||
|
||||
double precision,allocatable :: pt2(:,:)
|
||||
integer :: n_tasks, k, n_tasks_max
|
||||
integer :: n_tasks, k
|
||||
integer, allocatable :: i_generator(:), subset(:)
|
||||
|
||||
n_tasks_max = N_det_generators/100+1
|
||||
allocate(task_id(n_tasks_max), task(n_tasks_max))
|
||||
allocate(pt2(N_states,n_tasks_max), i_generator(n_tasks_max), subset(n_tasks_max))
|
||||
allocate(task_id(pt2_n_tasks_max), task(pt2_n_tasks_max))
|
||||
allocate(pt2(N_states,pt2_n_tasks_max), i_generator(pt2_n_tasks_max), subset(pt2_n_tasks_max))
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
|
||||
@ -47,7 +46,7 @@ subroutine run_pt2_slave(thread,iproc,energy)
|
||||
do while (.not.done)
|
||||
|
||||
n_tasks = max(1,n_tasks)
|
||||
n_tasks = min(n_tasks,n_tasks_max)
|
||||
n_tasks = min(n_tasks,pt2_n_tasks_max)
|
||||
|
||||
integer, external :: get_tasks_from_taskserver
|
||||
if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then
|
||||
@ -66,7 +65,7 @@ subroutine run_pt2_slave(thread,iproc,energy)
|
||||
do k=1,n_tasks
|
||||
pt2(:,k) = 0.d0
|
||||
buf%cur = 0
|
||||
call select_connected(i_generator(k),energy,pt2(1,k),buf,subset(k))
|
||||
call select_connected(i_generator(k),energy,pt2(1,k),buf,subset(k),pt2_F(i_generator(k)))
|
||||
enddo
|
||||
call wall_time(time1)
|
||||
|
||||
@ -78,7 +77,6 @@ subroutine run_pt2_slave(thread,iproc,energy)
|
||||
|
||||
! Try to adjust n_tasks around 1 second per job
|
||||
n_tasks = min(n_tasks,int( 1.d0*dble(n_tasks) / (time1 - time0 + 1.d-9)))+1
|
||||
! n_tasks = n_tasks+1
|
||||
end do
|
||||
|
||||
integer, external :: disconnect_from_taskserver
|
||||
@ -202,11 +200,4 @@ IRP_ENDIF
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, pt2_workload, (N_det_generators) ]
|
||||
integer :: i
|
||||
do i=1,N_det_generators
|
||||
pt2_workload(i) = dfloat(N_det_generators - i + 1)**2
|
||||
end do
|
||||
pt2_workload = pt2_workload / sum(pt2_workload)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
@ -1,122 +1,4 @@
|
||||
subroutine run_selection_slave(thread,iproc,energy)
|
||||
implicit none
|
||||
integer, intent(in) :: thread, iproc
|
||||
double precision, intent(in) :: energy(N_states_diag)
|
||||
call run_selection_slave_new(thread,iproc,energy)
|
||||
end
|
||||
|
||||
subroutine run_selection_slave_new(thread,iproc,energy)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: thread, iproc
|
||||
double precision, intent(in) :: energy(N_states_diag)
|
||||
integer :: rc, i, N
|
||||
logical :: buffer_ready
|
||||
|
||||
integer :: worker_id, ltask
|
||||
character*(512), allocatable :: task(:)
|
||||
integer, allocatable :: task_id(:)
|
||||
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
|
||||
integer(ZMQ_PTR), external :: new_zmq_push_socket
|
||||
integer(ZMQ_PTR) :: zmq_socket_push
|
||||
|
||||
type(selection_buffer) :: buf, buf2
|
||||
logical :: done
|
||||
|
||||
double precision,allocatable :: pt2(:,:)
|
||||
integer :: n_tasks, k, n_tasks_max
|
||||
integer, allocatable :: i_generator(:), subset(:)
|
||||
|
||||
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
|
||||
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
|
||||
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
|
||||
PROVIDE psi_bilinear_matrix_transp_order
|
||||
|
||||
buffer_ready = .False.
|
||||
n_tasks_max = N_det_generators/100+1
|
||||
allocate(task_id(n_tasks_max), task(n_tasks_max))
|
||||
allocate(pt2(N_states,n_tasks_max), i_generator(n_tasks_max), subset(n_tasks_max))
|
||||
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
|
||||
integer, external :: connect_to_taskserver
|
||||
if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
return
|
||||
endif
|
||||
|
||||
zmq_socket_push = new_zmq_push_socket(thread)
|
||||
|
||||
buf%N = 0
|
||||
n_tasks = 1
|
||||
call create_selection_buffer(0, 0, buf)
|
||||
done = .False.
|
||||
do while (.not.done)
|
||||
|
||||
n_tasks = max(1,n_tasks)
|
||||
n_tasks = min(n_tasks,n_tasks_max)
|
||||
|
||||
integer, external :: get_tasks_from_taskserver
|
||||
if (get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task, n_tasks) == -1) then
|
||||
exit
|
||||
endif
|
||||
done = task_id(n_tasks) == 0
|
||||
if (done) n_tasks = n_tasks-1
|
||||
if (n_tasks == 0) exit
|
||||
|
||||
do k=1,n_tasks
|
||||
read (task(k),*) subset(k), i_generator(k), N
|
||||
enddo
|
||||
|
||||
if(buf%N == 0) then
|
||||
! Only first time
|
||||
call create_selection_buffer(N, N*2, buf)
|
||||
call create_selection_buffer(N, N*2, buf2)
|
||||
buffer_ready = .True.
|
||||
endif
|
||||
|
||||
double precision :: time0, time1
|
||||
call wall_time(time0)
|
||||
do k=1,n_tasks
|
||||
pt2(:,k) = 0.d0
|
||||
buf%cur = 0
|
||||
call select_connected(i_generator(k),energy,pt2(1,k),buf,subset(k))
|
||||
enddo
|
||||
call wall_time(time1)
|
||||
|
||||
integer, external :: tasks_done_to_taskserver
|
||||
if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then
|
||||
done = .true.
|
||||
endif
|
||||
call sort_selection_buffer(buf)
|
||||
call merge_selection_buffers(buf,buf2)
|
||||
call push_selection_results(zmq_socket_push, pt2, buf, task_id, n_tasks)
|
||||
buf%mini = buf2%mini
|
||||
pt2(:,:) = 0d0
|
||||
buf%cur = 0
|
||||
|
||||
! ! Try to adjust n_tasks around 5 second per job
|
||||
! n_tasks = min(n_tasks,int( 5.d0 * dble(n_tasks) / (time1 - time0 + 1.d-9)))+1
|
||||
n_tasks = n_tasks+1
|
||||
end do
|
||||
|
||||
integer, external :: disconnect_from_taskserver
|
||||
if (disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id) == -1) then
|
||||
continue
|
||||
endif
|
||||
|
||||
call end_zmq_push_socket(zmq_socket_push,thread)
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call delete_selection_buffer(buf)
|
||||
|
||||
end
|
||||
|
||||
subroutine run_selection_slave_old(thread,iproc,energy)
|
||||
use f77_zmq
|
||||
use selection_types
|
||||
implicit none
|
||||
@ -177,7 +59,7 @@ subroutine run_selection_slave_old(thread,iproc,energy)
|
||||
else
|
||||
ASSERT (N == buf%N)
|
||||
end if
|
||||
call select_connected(i_generator,energy,pt2,buf,subset)
|
||||
call select_connected(i_generator,energy,pt2,buf,subset,pt2_F(i_generator))
|
||||
endif
|
||||
|
||||
integer, external :: task_done_to_taskserver
|
||||
|
||||
@ -1,14 +1,5 @@
|
||||
use bitmasks
|
||||
|
||||
BEGIN_PROVIDER [ integer, fragment_count ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of fragments for the deterministic part
|
||||
END_DOC
|
||||
fragment_count = (elec_alpha_num-n_core_orb)**2
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine assert(cond, msg)
|
||||
character(*), intent(in) :: msg
|
||||
logical, intent(in) :: cond
|
||||
@ -46,11 +37,11 @@ subroutine get_mask_phase(det, phasemask)
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine select_connected(i_generator,E0,pt2,b,subset)
|
||||
subroutine select_connected(i_generator,E0,pt2,b,subset,csubset)
|
||||
use bitmasks
|
||||
use selection_types
|
||||
implicit none
|
||||
integer, intent(in) :: i_generator, subset
|
||||
integer, intent(in) :: i_generator, subset, csubset
|
||||
type(selection_buffer), intent(inout) :: b
|
||||
double precision, intent(inout) :: pt2(N_states)
|
||||
integer :: k,l
|
||||
@ -71,7 +62,7 @@ subroutine select_connected(i_generator,E0,pt2,b,subset)
|
||||
particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) )
|
||||
particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) )
|
||||
enddo
|
||||
call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b,subset)
|
||||
call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b,subset,csubset)
|
||||
enddo
|
||||
deallocate(fock_diag_tmp)
|
||||
end subroutine
|
||||
@ -266,7 +257,7 @@ subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
|
||||
end
|
||||
|
||||
|
||||
subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf,subset)
|
||||
subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf,subset,csubset)
|
||||
use bitmasks
|
||||
use selection_types
|
||||
implicit none
|
||||
@ -274,7 +265,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
|
||||
! WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: i_generator, subset
|
||||
integer, intent(in) :: i_generator, subset, csubset
|
||||
integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2)
|
||||
double precision, intent(in) :: fock_diag_tmp(mo_tot_num)
|
||||
double precision, intent(in) :: E0(N_states)
|
||||
@ -298,8 +289,6 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
|
||||
integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
|
||||
allocate (preinteresting_det(N_int,2,N_det))
|
||||
|
||||
PROVIDE fragment_count
|
||||
|
||||
monoAdo = .true.
|
||||
monoBdo = .true.
|
||||
|
||||
@ -571,7 +560,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
|
||||
end if
|
||||
|
||||
maskInd += 1
|
||||
if(subset == 0 .or. mod(maskInd, fragment_count) == (subset-1)) then
|
||||
if(mod(maskInd, csubset) == (subset-1)) then
|
||||
|
||||
call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
|
||||
if(fullMatch) cycle
|
||||
|
||||
@ -14,7 +14,7 @@ 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 n_states_diag
|
||||
PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count ci_energy mpi_master zmq_state zmq_context
|
||||
PROVIDE pt2_e0_denominator mo_tot_num N_int ci_energy mpi_master zmq_state zmq_context
|
||||
PROVIDE psi_det psi_coef
|
||||
end
|
||||
|
||||
|
||||
@ -14,7 +14,7 @@ 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 n_states_diag
|
||||
PROVIDE pt2_e0_denominator mo_tot_num N_int fragment_count ci_energy mpi_master zmq_state zmq_context
|
||||
PROVIDE pt2_e0_denominator mo_tot_num N_int ci_energy mpi_master zmq_state zmq_context
|
||||
PROVIDE psi_det psi_coef
|
||||
end
|
||||
|
||||
|
||||
@ -12,15 +12,13 @@ subroutine ZMQ_selection(N_in, pt2)
|
||||
double precision, intent(out) :: pt2(N_states)
|
||||
|
||||
|
||||
PROVIDE fragment_count
|
||||
|
||||
N = max(N_in,1)
|
||||
if (.True.) then
|
||||
PROVIDE pt2_e0_denominator nproc
|
||||
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
|
||||
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
|
||||
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
|
||||
PROVIDE psi_bilinear_matrix_transp_order fragment_count
|
||||
PROVIDE psi_bilinear_matrix_transp_order
|
||||
|
||||
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')
|
||||
|
||||
@ -60,9 +58,8 @@ subroutine ZMQ_selection(N_in, pt2)
|
||||
task = ' '
|
||||
|
||||
do i= 1, N_det_generators
|
||||
! /!\ Fragments don't work
|
||||
! if (i>-ishft(N_det_generators,-2)) then
|
||||
write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') 0, i, N
|
||||
do j=1,pt2_F(pt2_J(i))
|
||||
write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, pt2_J(i), N
|
||||
ipos += 30
|
||||
if (ipos > 63970) then
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
|
||||
@ -70,18 +67,7 @@ subroutine ZMQ_selection(N_in, pt2)
|
||||
endif
|
||||
ipos=1
|
||||
endif
|
||||
! else
|
||||
! do j=1,fragment_count
|
||||
! write(task(ipos:ipos+30),'(I9,1X,I9,1X,I9,''|'')') j, i, N
|
||||
! ipos += 30
|
||||
! if (ipos > 63970) then
|
||||
! if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
|
||||
! stop 'Unable to add task to task server'
|
||||
! endif
|
||||
! ipos=1
|
||||
! endif
|
||||
! end do
|
||||
! endif
|
||||
end do
|
||||
enddo
|
||||
if (ipos > 1) then
|
||||
if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then
|
||||
|
||||
@ -13,7 +13,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
|
||||
N_det_generators = N_det
|
||||
do i=1,N_det
|
||||
norm = norm + psi_average_norm_contrib_sorted(i)
|
||||
if (norm >= threshold_generators) then
|
||||
if (norm > threshold_generators+1d-10) then
|
||||
N_det_generators = i
|
||||
exit
|
||||
endif
|
||||
@ -29,7 +29,6 @@ END_PROVIDER
|
||||
! For Single reference wave functions, the generator is the
|
||||
! Hartree-Fock determinant
|
||||
END_DOC
|
||||
integer :: i, k
|
||||
psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
|
||||
psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states)
|
||||
|
||||
@ -44,7 +43,6 @@ END_PROVIDER
|
||||
! For Single reference wave functions, the generator is the
|
||||
! Hartree-Fock determinant
|
||||
END_DOC
|
||||
integer :: i, k
|
||||
psi_det_sorted_gen = psi_det_sorted
|
||||
psi_coef_sorted_gen = psi_coef_sorted
|
||||
psi_det_sorted_gen_order = psi_det_sorted_order
|
||||
|
||||
@ -17,12 +17,6 @@ doc: Stop stochastic PT2 when the relative error is smaller than PT2_relative_er
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.001
|
||||
|
||||
[PT2_absolute_error]
|
||||
type: Threshold
|
||||
doc: Stop stochastic PT2 when the statistical error is smaller than PT2_absolute_error
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.00001
|
||||
|
||||
[correlation_energy_ratio_max]
|
||||
type: Normalized_float
|
||||
doc: The selection process stops at a fixed correlation ratio (useful for getting same accuracy between molecules)
|
||||
|
||||
@ -10,3 +10,9 @@ doc: Maximum number of dressed CI iterations
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 10
|
||||
|
||||
[dress_relative_error]
|
||||
type: Normalized_float
|
||||
doc: Stop stochastic PT2 when the relative error is smaller than PT2_relative_error
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.001
|
||||
|
||||
|
||||
@ -2,10 +2,10 @@ use bitmasks
|
||||
|
||||
|
||||
|
||||
subroutine alpha_callback(delta_ij_loc, i_generator, subset,iproc)
|
||||
subroutine alpha_callback(delta_ij_loc, i_generator, subset, csubset, iproc)
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer, intent(in) :: i_generator, subset
|
||||
integer, intent(in) :: i_generator, subset, csubset
|
||||
double precision,intent(inout) :: delta_ij_loc(N_states,N_det,2)
|
||||
integer, intent(in) :: iproc
|
||||
|
||||
@ -15,7 +15,7 @@ subroutine alpha_callback(delta_ij_loc, i_generator, subset,iproc)
|
||||
|
||||
|
||||
do l=1,N_generators_bitmask
|
||||
call generate_singles_and_doubles(delta_ij_loc, i_generator,l,subset,iproc)
|
||||
call generate_singles_and_doubles(delta_ij_loc,i_generator,l,subset,csubset,iproc)
|
||||
enddo
|
||||
end subroutine
|
||||
|
||||
@ -34,7 +34,7 @@ BEGIN_PROVIDER [ integer, psi_from_sorted_gen, (N_det) ]
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index, subset, iproc)
|
||||
subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index, subset, csubset, iproc)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -42,7 +42,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
|
||||
END_DOC
|
||||
|
||||
double precision,intent(inout) :: delta_ij_loc(N_states,N_det,2)
|
||||
integer, intent(in) :: i_generator, subset, bitmask_index
|
||||
integer, intent(in) :: i_generator, subset, csubset, bitmask_index
|
||||
integer, intent(in) :: iproc
|
||||
|
||||
|
||||
@ -66,11 +66,11 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
|
||||
integer(bit_kind), allocatable:: preinteresting_det(:,:,:)
|
||||
integer ,allocatable :: abuf(:), labuf(:)
|
||||
|
||||
allocate(abuf(0:N_det*6), labuf(0:N_det))
|
||||
allocate(abuf(N_det*6), labuf(N_det))
|
||||
allocate(preinteresting_det(N_int,2,N_det))
|
||||
|
||||
PROVIDE fragment_count
|
||||
|
||||
maskInd = -1
|
||||
|
||||
monoAdo = .true.
|
||||
monoBdo = .true.
|
||||
@ -193,7 +193,6 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
|
||||
allocate(counted(mo_tot_num, mo_tot_num), countedOrb(mo_tot_num, 2))
|
||||
allocate (indexes(0:mo_tot_num, 0:mo_tot_num))
|
||||
allocate (indexes_end(0:mo_tot_num, 0:mo_tot_num))
|
||||
maskInd = -1
|
||||
integer :: nb_count
|
||||
do s1=1,2
|
||||
do i1=N_holes(s1),1,-1 ! Generate low excitations first
|
||||
@ -345,7 +344,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
|
||||
end if
|
||||
|
||||
maskInd += 1
|
||||
if(subset == 0 .or. mod(maskInd, fragment_count) == (subset-1)) then
|
||||
if(mod(maskInd, csubset) == (subset-1)) then
|
||||
|
||||
call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
|
||||
if(fullMatch) cycle
|
||||
@ -387,7 +386,7 @@ subroutine alpha_callback_mask(delta_ij_loc, i_gen, sp, mask, bannedOrb, banned,
|
||||
integer(bit_kind), allocatable :: det_minilist(:,:,:)
|
||||
|
||||
|
||||
allocate(abuf(0:siz), labuf(0:N_det), putten(N_det), det_minilist(N_int, 2, N_det))
|
||||
allocate(abuf(siz), labuf(N_det), putten(N_det), det_minilist(N_int, 2, N_det))
|
||||
|
||||
do i=1,siz
|
||||
abuf(i) = psi_from_sorted_gen(rabuf(i))
|
||||
@ -638,7 +637,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, indexes, ab
|
||||
integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel)
|
||||
logical, intent(inout) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num, 2)
|
||||
integer, intent(inout) :: indexes(0:mo_tot_num, 0:mo_tot_num)
|
||||
integer, intent(inout) :: abuf(0:*)
|
||||
integer, intent(inout) :: abuf(*)
|
||||
integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt, s
|
||||
integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
|
||||
integer :: phasemask(2,N_int*bit_kind_size)
|
||||
@ -704,7 +703,7 @@ subroutine get_d2(i_gen, gen, banned, bannedOrb, indexes, abuf, mask, h, p, sp)
|
||||
implicit none
|
||||
|
||||
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
|
||||
integer, intent(inout) :: abuf(0:*)
|
||||
integer, intent(inout) :: abuf(*)
|
||||
integer, intent(in) :: i_gen
|
||||
logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
|
||||
integer, intent(inout) :: indexes(0:mo_tot_num, 0:mo_tot_num)
|
||||
@ -832,7 +831,7 @@ subroutine get_d1(i_gen, gen, banned, bannedOrb, indexes, abuf, mask, h, p, sp)
|
||||
implicit none
|
||||
|
||||
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
|
||||
integer, intent(inout) :: abuf(0:*)
|
||||
integer, intent(inout) :: abuf(*)
|
||||
integer,intent(in) :: i_gen
|
||||
logical, intent(in) :: bannedOrb(mo_tot_num, 2), banned(mo_tot_num, mo_tot_num,2)
|
||||
integer(bit_kind) :: det(N_int, 2)
|
||||
|
||||
@ -29,8 +29,6 @@ subroutine run_dressing(N_st,energy)
|
||||
delta_E = 1.d0
|
||||
iteration = 0
|
||||
do iteration=1,n_it_dress_max
|
||||
N_det_delta_ij = N_det
|
||||
touch N_det_delta_ij
|
||||
print *, '==============================================='
|
||||
print *, 'Iteration', iteration, '/', n_it_dress_max
|
||||
print *, '==============================================='
|
||||
@ -40,13 +38,11 @@ subroutine run_dressing(N_st,energy)
|
||||
do i=1,N_st
|
||||
print *, i, psi_energy(i)+nuclear_repulsion
|
||||
enddo
|
||||
!print *, "DELTA IJ", delta_ij(1,1,1)
|
||||
PROVIDE delta_ij_tmp
|
||||
if(.true.) call delta_ij_done()
|
||||
print *, 'Dressed energy <Psi|H+Delta|Psi>'
|
||||
do i=1,N_st
|
||||
print *, i, ci_energy_dressed(i)
|
||||
enddo
|
||||
energy(1:N_st) = ci_energy_dressed(1:N_st)
|
||||
call diagonalize_ci_dressed
|
||||
E_new = sum(psi_energy(:))
|
||||
|
||||
@ -56,7 +52,6 @@ subroutine run_dressing(N_st,energy)
|
||||
call write_double(6,delta_E,"delta_E (undressed)")
|
||||
delta_E = dabs(delta_E)
|
||||
call save_wavefunction
|
||||
! call ezfio_set_dress_zmq_energy(ci_energy_dressed(1))
|
||||
if (delta_E < thresh_dress) then
|
||||
exit
|
||||
endif
|
||||
@ -67,10 +62,9 @@ subroutine run_dressing(N_st,energy)
|
||||
enddo
|
||||
print *, 'Dressed energy <Psi|H+Delta|Psi>'
|
||||
do i=1,N_st
|
||||
print *, i, ci_energy_dressed(i)+nuclear_repulsion
|
||||
print *, i, ci_energy_dressed(i)
|
||||
enddo
|
||||
endif
|
||||
|
||||
if(.true.) energy(1:N_st) = 0d0 ! ci_energy_dressed(1:N_st)
|
||||
end
|
||||
|
||||
|
||||
@ -50,9 +50,7 @@ subroutine run_wf
|
||||
else if (zmq_state(:5) == 'dress') then
|
||||
! Dress
|
||||
! ---------
|
||||
!call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states)
|
||||
if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
|
||||
!TOUCH psi_det
|
||||
if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle
|
||||
if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle
|
||||
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -1,71 +1,9 @@
|
||||
use bitmasks
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, N_dress_teeth ]
|
||||
N_dress_teeth = 10
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, dress_norm_acc, (0:N_det, N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, dress_norm, (0:N_det, N_states) ]
|
||||
&BEGIN_PROVIDER [ double precision, dress_teeth_size, (0:N_det, N_states) ]
|
||||
&BEGIN_PROVIDER [ integer, dress_teeth, (0:N_dress_teeth+1, N_states) ]
|
||||
implicit none
|
||||
integer :: i, j, st, nt
|
||||
double precision :: norm_sto, jump, norm_mwen, norm_loc
|
||||
|
||||
if(N_states /= 1) stop "dress_sto may not work with N_states /= 1"
|
||||
|
||||
do st=1,N_states
|
||||
dress_teeth(0,st) = 1
|
||||
norm_sto = 1d0
|
||||
do i=1,N_det
|
||||
dress_teeth(1,st) = i
|
||||
jump = (1d0 / dfloat(N_dress_teeth)) * norm_sto
|
||||
if(psi_coef_generators(i,1)**2 < jump / 2d0) exit
|
||||
norm_sto -= psi_coef_generators(i,1)**2
|
||||
end do
|
||||
|
||||
norm_loc = 0d0
|
||||
dress_norm_acc(0,st) = 0d0
|
||||
nt = 1
|
||||
|
||||
do i=1,dress_teeth(1,st)-1
|
||||
dress_norm_acc(i,st) = dress_norm_acc(i-1,st) + psi_coef_generators(i,st)**2
|
||||
end do
|
||||
|
||||
do i=dress_teeth(1,st), N_det_generators!-dress_teeth(1,st)+1
|
||||
norm_mwen = psi_coef_generators(i,st)**2!-1+dress_teeth(1,st),st)**2
|
||||
dress_norm_acc(i,st) = dress_norm_acc(i-1,st) + norm_mwen
|
||||
norm_loc += norm_mwen
|
||||
if(norm_loc > (jump*dfloat(nt))) then
|
||||
nt = nt + 1
|
||||
dress_teeth(nt,st) = i
|
||||
end if
|
||||
end do
|
||||
if(nt > N_dress_teeth+1) then
|
||||
print *, "foireouse dress_teeth", nt, dress_teeth(nt,st), N_det
|
||||
stop
|
||||
end if
|
||||
|
||||
dress_teeth(N_dress_teeth+1,st) = N_det+1
|
||||
norm_loc = 0d0
|
||||
do i=N_dress_teeth, 0, -1
|
||||
dress_teeth_size(i,st) = dress_norm_acc(dress_teeth(i+1,st)-1,st) - dress_norm_acc(dress_teeth(i,st)-1, st)
|
||||
dress_norm_acc(dress_teeth(i,st):dress_teeth(i+1,st)-1,st) -= dress_norm_acc(dress_teeth(i,st)-1, st)
|
||||
dress_norm_acc(dress_teeth(i,st):dress_teeth(i+1,st)-1,st) = &
|
||||
dress_norm_acc(dress_teeth(i,st):dress_teeth(i+1,st)-1,st) / dress_teeth_size(i,st)
|
||||
dress_norm(dress_teeth(i,st), st) = dress_norm_acc(dress_teeth(i,st), st)
|
||||
do j=dress_teeth(i,st)+1, dress_teeth(i+1,1)-1
|
||||
dress_norm(j,1) = dress_norm_acc(j, st) - dress_norm_acc(j-1, st)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer , N_det_delta_ij ]
|
||||
implicit none
|
||||
N_det_delta_ij = 1
|
||||
N_det_delta_ij = N_det
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_ij, (N_states, N_det, 2) ]
|
||||
@ -83,36 +21,23 @@ BEGIN_PROVIDER [ double precision, delta_ij_tmp, (N_states,N_det_delta_ij,2) ]
|
||||
integer :: i,j,k
|
||||
|
||||
double precision, allocatable :: dress(:), del(:,:), del_s2(:,:)
|
||||
double precision :: E_CI_before(N_states), relative_error
|
||||
|
||||
! prevents re-providing if delta_ij_tmp is
|
||||
! just being copied
|
||||
if(N_det_delta_ij == N_det) then
|
||||
double precision :: E_CI_before(N_states)
|
||||
integer :: cnt = 0
|
||||
|
||||
allocate(dress(N_states), del(N_states, N_det_delta_ij), del_s2(N_states, N_det_delta_ij))
|
||||
|
||||
delta_ij_tmp = 0d0
|
||||
|
||||
E_CI_before(:) = psi_energy(:) + nuclear_repulsion
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1.d0
|
||||
SOFT_TOUCH threshold_selectors threshold_generators
|
||||
! if(errr /= 0d0) then
|
||||
! errr = errr / 2d0
|
||||
! else
|
||||
! errr = 1d-4
|
||||
! end if
|
||||
relative_error = 1.d-3
|
||||
E_CI_before(:) = dress_E0_denominator(:) + nuclear_repulsion
|
||||
|
||||
call write_double(6,relative_error,"Relative error for the stochastic algorithm")
|
||||
call write_double(6,dress_relative_error,"Convergence of the stochastic algorithm")
|
||||
|
||||
call ZMQ_dress(E_CI_before, dress, del, del_s2, abs(relative_error), N_det_delta_ij)
|
||||
call ZMQ_dress(E_CI_before, dress, del, del_s2, abs(dress_relative_error), N_det_delta_ij)
|
||||
delta_ij_tmp(:,:,1) = del(:,:)
|
||||
delta_ij_tmp(:,:,2) = del_s2(:,:)
|
||||
|
||||
|
||||
deallocate(dress, del, del_s2)
|
||||
end if
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
@ -1,13 +1,5 @@
|
||||
use bitmasks
|
||||
|
||||
BEGIN_PROVIDER [ integer, fragment_count ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of fragments for the deterministic part
|
||||
END_DOC
|
||||
fragment_count = 1
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine run_dress_slave(thread,iproce,energy)
|
||||
use f77_zmq
|
||||
@ -16,10 +8,10 @@ subroutine run_dress_slave(thread,iproce,energy)
|
||||
|
||||
double precision, intent(in) :: energy(N_states_diag)
|
||||
integer, intent(in) :: thread, iproce
|
||||
integer :: rc, i, subset, i_generator
|
||||
integer :: rc, i, j, subset, i_generator
|
||||
|
||||
integer :: worker_id, task_id, ctask, ltask
|
||||
character*(5120) :: task
|
||||
integer :: worker_id, ctask, ltask
|
||||
character*(512) :: task(Nproc)
|
||||
|
||||
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
|
||||
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
|
||||
@ -27,70 +19,64 @@ subroutine run_dress_slave(thread,iproce,energy)
|
||||
integer(ZMQ_PTR), external :: new_zmq_push_socket
|
||||
integer(ZMQ_PTR) :: zmq_socket_push
|
||||
|
||||
logical :: done
|
||||
|
||||
double precision,allocatable :: dress_detail(:)
|
||||
integer :: ind
|
||||
|
||||
double precision,allocatable :: delta_ij_loc(:,:,:)
|
||||
integer :: h,p,n,i_state
|
||||
logical :: ok
|
||||
|
||||
integer, allocatable :: int_buf(:)
|
||||
double precision, allocatable :: double_buf(:)
|
||||
integer(bit_kind), allocatable :: det_buf(:,:,:)
|
||||
integer :: N_buf(3)
|
||||
logical :: last
|
||||
double precision,allocatable :: breve_delta_m(:,:,:)
|
||||
integer :: i_state,m,l,t,p,sum_f
|
||||
!integer, external :: omp_get_thread_num
|
||||
double precision, allocatable :: delta_det(:,:,:,:), cp(:,:,:,:)
|
||||
integer :: toothMwen
|
||||
logical :: fracted
|
||||
double precision, allocatable :: delta_det(:,:,:,:), cp(:,:,:,:), edI(:)
|
||||
double precision, allocatable :: edI_task(:)
|
||||
integer, allocatable :: edI_index(:), edI_taskID(:)
|
||||
integer :: n_tasks
|
||||
|
||||
integer :: iproc
|
||||
integer, allocatable :: f(:)
|
||||
integer :: cp_sent, cp_done
|
||||
integer :: cp_max(Nproc)
|
||||
integer :: will_send, task_id, purge_task_id, ntask_buf
|
||||
integer, allocatable :: task_buf(:)
|
||||
integer(kind=OMP_LOCK_KIND) :: lck_det(0:pt2_N_teeth+1)
|
||||
integer(kind=OMP_LOCK_KIND) :: lck_sto(0:dress_N_cp+1), sending, getting_task
|
||||
double precision :: fac
|
||||
|
||||
|
||||
double precision :: ending(1)
|
||||
integer, external :: zmq_get_dvector
|
||||
! double precision, external :: omp_get_wtime
|
||||
double precision :: time, time0
|
||||
integer :: ntask_tbd, task_tbd(Nproc), i_gen_tbd(Nproc), subset_tbd(Nproc)
|
||||
! if(iproce /= 0) stop "RUN DRESS SLAVE is OMP"
|
||||
|
||||
allocate(delta_det(N_states, N_det, 0:comb_teeth+1, 2))
|
||||
allocate(cp(N_states, N_det, N_cp, 2))
|
||||
delta_det = 0d9
|
||||
allocate(delta_det(N_states, N_det, 0:pt2_N_teeth+1, 2))
|
||||
allocate(cp(N_states, N_det, dress_N_cp, 2))
|
||||
allocate(edI(N_det_generators), f(N_det_generators))
|
||||
allocate(edI_index(N_det_generators), edI_task(N_det_generators))
|
||||
edI = 0d0
|
||||
f = 0
|
||||
delta_det = 0d0
|
||||
cp = 0d0
|
||||
task = CHAR(0)
|
||||
|
||||
|
||||
task(:) = CHAR(0)
|
||||
|
||||
|
||||
|
||||
integer :: iproc, cur_cp, done_for(0:N_cp)
|
||||
integer, allocatable :: tasks(:)
|
||||
integer :: lastCp(Nproc)
|
||||
integer :: lastSent, lastSendable
|
||||
logical :: send
|
||||
integer(kind=OMP_LOCK_KIND) :: lck_det(0:comb_teeth+1)
|
||||
integer(kind=OMP_LOCK_KIND) :: lck_sto(0:N_cp+1)
|
||||
|
||||
do i=0,N_cp+1
|
||||
call omp_init_lock(sending)
|
||||
call omp_init_lock(getting_task)
|
||||
do i=0,dress_N_cp+1
|
||||
call omp_init_lock(lck_sto(i))
|
||||
end do
|
||||
do i=0,comb_teeth+1
|
||||
do i=0,pt2_N_teeth+1
|
||||
call omp_init_lock(lck_det(i))
|
||||
end do
|
||||
|
||||
lastCp = 0
|
||||
lastSent = 0
|
||||
send = .false.
|
||||
done_for = 0
|
||||
|
||||
double precision :: hij, sij
|
||||
!call i_h_j_s2(psi_det(1,1,1),psi_det(1,1,2),N_int,hij, sij)
|
||||
|
||||
hij = dress_E0_denominator(1) !PROVIDE BEFORE OMP PARALLEL
|
||||
cp_done = 0
|
||||
cp_sent = 0
|
||||
will_send = 0
|
||||
|
||||
double precision :: hij, sij, tmp
|
||||
purge_task_id = 0
|
||||
provide psi_energy
|
||||
ending(1) = dble(dress_N_cp+1)
|
||||
ntask_tbd = 0
|
||||
!$OMP PARALLEL DEFAULT(SHARED) &
|
||||
!$OMP PRIVATE(int_buf, double_buf, det_buf, delta_ij_loc, task, task_id) &
|
||||
!$OMP PRIVATE(lastSendable, toothMwen, fracted, fac) &
|
||||
!$OMP PRIVATE(i, cur_cp, send, i_generator, subset, iproc, N_buf) &
|
||||
!$OMP PRIVATE(zmq_to_qp_run_socket, zmq_socket_push, worker_id)
|
||||
|
||||
!$OMP PRIVATE(breve_delta_m, task_id) &
|
||||
!$OMP PRIVATE(tmp,fac,m,l,t,sum_f,n_tasks) &
|
||||
!$OMP PRIVATE(i,p,will_send, i_generator, subset, iproc) &
|
||||
!$OMP PRIVATE(zmq_to_qp_run_socket, zmq_socket_push, worker_id) &
|
||||
!$OMP PRIVATE(task_buf, ntask_buf,time, time0)
|
||||
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
|
||||
zmq_socket_push = new_zmq_push_socket(thread)
|
||||
call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread)
|
||||
@ -99,273 +85,244 @@ subroutine run_dress_slave(thread,iproce,energy)
|
||||
call end_zmq_push_socket(zmq_socket_push,thread)
|
||||
stop "WORKER -1"
|
||||
end if
|
||||
|
||||
|
||||
iproc = omp_get_thread_num()+1
|
||||
allocate(int_buf(N_dress_int_buffer))
|
||||
allocate(double_buf(N_dress_double_buffer))
|
||||
allocate(det_buf(N_int, 2, N_dress_det_buffer))
|
||||
allocate(delta_ij_loc(N_states,N_det,2))
|
||||
do
|
||||
call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task)
|
||||
task = task//" 0"
|
||||
if(task_id == 0) exit
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
allocate(breve_delta_m(N_states,N_det,2))
|
||||
allocate(task_buf(pt2_n_tasks_max))
|
||||
ntask_buf = 0
|
||||
|
||||
if(iproc==1) then
|
||||
call push_dress_results(zmq_socket_push, 0, 0, edI_task, edI_index, breve_delta_m, task_buf, ntask_buf)
|
||||
end if
|
||||
|
||||
do while(cp_done > cp_sent .or. m /= dress_N_cp+1)
|
||||
call omp_set_lock(getting_task)
|
||||
if(ntask_tbd == 0) then
|
||||
ntask_tbd = size(task_tbd)
|
||||
call get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id, task_tbd, task, ntask_tbd)
|
||||
!task = task//" 0"
|
||||
end if
|
||||
|
||||
task_id = task_tbd(1)
|
||||
if(task_id /= 0) then
|
||||
read (task,*) subset, i_generator
|
||||
|
||||
!$OMP ATOMIC
|
||||
done_for(done_cp_at_det(i_generator)) += 1
|
||||
! print *, "IGEN", i_generator, done_cp_at_det(i_generator)
|
||||
delta_ij_loc(:,:,:) = 0d0
|
||||
call generator_start(i_generator, iproc)
|
||||
call alpha_callback(delta_ij_loc, i_generator, subset, iproc)
|
||||
call generator_done(i_generator, int_buf, double_buf, det_buf, N_buf, iproc)
|
||||
|
||||
do i=1,N_cp
|
||||
fac = cps(i_generator, i) * dress_weight_inv(i_generator) * comb_step
|
||||
if(fac == 0d0) cycle
|
||||
call omp_set_lock(lck_sto(i))
|
||||
cp(:,:,i,1) += (delta_ij_loc(:,:,1) * fac)
|
||||
cp(:,:,i,2) += (delta_ij_loc(:,:,2) * fac)
|
||||
call omp_unset_lock(lck_sto(i))
|
||||
read (task(1),*) subset, i_generator
|
||||
do i=1,size(task_tbd)-1
|
||||
task_tbd(i) = task_tbd(i+1)
|
||||
task(i) = task(i+1)
|
||||
end do
|
||||
|
||||
|
||||
toothMwen = tooth_of_det(i_generator)
|
||||
fracted = (toothMwen /= 0)
|
||||
if(fracted) fracted = (i_generator == first_det_of_teeth(toothMwen))
|
||||
if(fracted) then
|
||||
call omp_set_lock(lck_det(toothMwen))
|
||||
call omp_set_lock(lck_det(toothMwen-1))
|
||||
delta_det(:,:,toothMwen-1, 1) += delta_ij_loc(:,:,1) * (1d0-fractage(toothMwen))
|
||||
delta_det(:,:,toothMwen-1, 2) += delta_ij_loc(:,:,2) * (1d0-fractage(toothMwen))
|
||||
delta_det(:,:,toothMwen , 1) += delta_ij_loc(:,:,1) * (fractage(toothMwen))
|
||||
delta_det(:,:,toothMwen , 2) += delta_ij_loc(:,:,2) * (fractage(toothMwen))
|
||||
call omp_unset_lock(lck_det(toothMwen))
|
||||
call omp_unset_lock(lck_det(toothMwen-1))
|
||||
m = dress_P(i_generator)
|
||||
ntask_tbd -= 1
|
||||
else
|
||||
call omp_set_lock(lck_det(toothMwen))
|
||||
delta_det(:,:,toothMwen , 1) += delta_ij_loc(:,:,1)
|
||||
delta_det(:,:,toothMwen , 2) += delta_ij_loc(:,:,2)
|
||||
call omp_unset_lock(lck_det(toothMwen))
|
||||
end if
|
||||
call push_dress_results(zmq_socket_push, i_generator, -1, delta_ij_loc, int_buf, double_buf, det_buf, N_buf, task_id)
|
||||
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
|
||||
lastCp(iproc) = done_cp_at_det(i_generator)
|
||||
m = dress_N_cp + 1
|
||||
i= zmq_get_dvector(zmq_to_qp_run_socket, worker_id, "ending", ending, 1)
|
||||
end if
|
||||
call omp_unset_lock(getting_task)
|
||||
will_send = 0
|
||||
|
||||
!$OMP CRITICAL
|
||||
send = .false.
|
||||
lastSendable = N_cp*2
|
||||
do i=1,Nproc
|
||||
lastSendable = min(lastCp(i), lastSendable)
|
||||
end do
|
||||
lastSendable -= 1
|
||||
if(lastSendable > lastSent .or. (lastSendable == N_cp-1 .and. lastSent /= N_cp-1)) then
|
||||
lastSent = lastSendable
|
||||
cur_cp = lastSent
|
||||
send = .true.
|
||||
cp_max(iproc) = m
|
||||
cp_done = minval(cp_max)-1
|
||||
if(cp_done > cp_sent) then
|
||||
will_send = cp_sent + 1
|
||||
cp_sent = will_send
|
||||
end if
|
||||
if(purge_task_id == 0) then
|
||||
purge_task_id = task_id
|
||||
task_id = 0
|
||||
else if(task_id /= 0) then
|
||||
ntask_buf += 1
|
||||
task_buf(ntask_buf) = task_id
|
||||
end if
|
||||
!$OMP END CRITICAL
|
||||
|
||||
if(send) then
|
||||
N_buf = (/0,1,0/)
|
||||
|
||||
delta_ij_loc = 0d0
|
||||
if(cur_cp < 1) stop "cur_cp < 1"
|
||||
do i=1,cur_cp
|
||||
delta_ij_loc(:,:,1) += cp(:,:,i,1)
|
||||
delta_ij_loc(:,:,2) += cp(:,:,i,2)
|
||||
if(will_send /= 0 .and. will_send <= int(ending(1))) then
|
||||
call omp_set_lock(sending)
|
||||
n_tasks = 0
|
||||
sum_f = 0
|
||||
do i=1,N_det_generators
|
||||
if(dress_P(i) <= will_send) sum_f = sum_f + f(i)
|
||||
if(dress_P(i) == will_send .and. f(i) /= 0) then
|
||||
n_tasks += 1
|
||||
edI_task(n_tasks) = edI(i)
|
||||
edI_index(n_tasks) = i
|
||||
end if
|
||||
end do
|
||||
|
||||
delta_ij_loc(:,:,:) = delta_ij_loc(:,:,:) / cps_N(cur_cp)
|
||||
do i=cp_first_tooth(cur_cp)-1,0,-1
|
||||
delta_ij_loc(:,:,1) = delta_ij_loc(:,:,1) +delta_det(:,:,i,1)
|
||||
delta_ij_loc(:,:,2) = delta_ij_loc(:,:,2) +delta_det(:,:,i,2)
|
||||
end do
|
||||
call push_dress_results(zmq_socket_push, done_for(cur_cp), cur_cp, delta_ij_loc, int_buf, double_buf, det_buf, N_buf, -1)
|
||||
call push_dress_results(zmq_socket_push, will_send, sum_f, edI_task, edI_index, breve_delta_m, 0, n_tasks)
|
||||
call omp_unset_lock(sending)
|
||||
end if
|
||||
|
||||
if(task_id == 0) exit
|
||||
if(m /= dress_N_cp+1) then
|
||||
!UPDATE i_generator
|
||||
|
||||
breve_delta_m(:,:,:) = 0d0
|
||||
call generator_start(i_generator, iproc)
|
||||
time0 = omp_get_wtime()
|
||||
call alpha_callback(breve_delta_m, i_generator, subset, pt2_F(i_generator), iproc)
|
||||
time = omp_get_wtime()
|
||||
!print '(I0.11, I4, A12, F12.3)', i_generator, subset, "GREPMETIME", time-time0
|
||||
t = dress_T(i_generator)
|
||||
|
||||
call omp_set_lock(lck_det(t))
|
||||
do j=1,N_det
|
||||
do i=1,N_states
|
||||
delta_det(i,j,t, 1) = delta_det(i,j,t, 1) + breve_delta_m(i,j,1)
|
||||
delta_det(i,j,t, 2) = delta_det(i,j,t, 2) + breve_delta_m(i,j,2)
|
||||
enddo
|
||||
enddo
|
||||
call omp_unset_lock(lck_det(t))
|
||||
|
||||
do p=1,dress_N_cp
|
||||
if(dress_e(i_generator, p) /= 0d0) then
|
||||
fac = dress_e(i_generator, p)
|
||||
call omp_set_lock(lck_sto(p))
|
||||
do j=1,N_det
|
||||
do i=1,N_states
|
||||
cp(i,j,p,1) = cp(i,j,p,1) + breve_delta_m(i,j,1) * fac
|
||||
cp(i,j,p,2) = cp(i,j,p,2) + breve_delta_m(i,j,2) * fac
|
||||
enddo
|
||||
enddo
|
||||
call omp_unset_lock(lck_sto(p))
|
||||
end if
|
||||
end do
|
||||
|
||||
tmp = 0d0
|
||||
do i=N_det,1,-1
|
||||
tmp += psi_coef(i, dress_stoch_istate)*breve_delta_m(dress_stoch_istate, i, 1)
|
||||
end do
|
||||
!$OMP ATOMIC
|
||||
edI(i_generator) += tmp
|
||||
!$OMP ATOMIC
|
||||
f(i_generator) += 1
|
||||
!push bidon
|
||||
if(ntask_buf == size(task_buf)) then
|
||||
call push_dress_results(zmq_socket_push, 0, 0, edI_task, edI_index, breve_delta_m, task_buf, ntask_buf)
|
||||
ntask_buf = 0
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
!$OMP BARRIER
|
||||
if(ntask_buf /= 0) then
|
||||
call push_dress_results(zmq_socket_push, 0, 0, edI_task, edI_index, breve_delta_m, task_buf, ntask_buf)
|
||||
ntask_buf = 0
|
||||
end if
|
||||
!$OMP SINGLE
|
||||
if(purge_task_id /= 0) then
|
||||
do while(int(ending(1)) == dress_N_cp+1)
|
||||
call sleep(1)
|
||||
i= zmq_get_dvector(zmq_to_qp_run_socket, worker_id, "ending", ending, 1)
|
||||
end do
|
||||
|
||||
will_send = int(ending(1))
|
||||
breve_delta_m = 0d0
|
||||
|
||||
do l=will_send, 1,-1
|
||||
breve_delta_m(:,:,1) += cp(:,:,l,1)
|
||||
breve_delta_m(:,:,2) += cp(:,:,l,2)
|
||||
end do
|
||||
|
||||
breve_delta_m(:,:,:) = breve_delta_m(:,:,:) / dress_M_m(will_send)
|
||||
|
||||
do t=dress_dot_t(will_send)-1,0,-1
|
||||
breve_delta_m(:,:,1) = breve_delta_m(:,:,1) + delta_det(:,:,t,1)
|
||||
breve_delta_m(:,:,2) = breve_delta_m(:,:,2) + delta_det(:,:,t,2)
|
||||
end do
|
||||
|
||||
sum_f = 0
|
||||
do i=1,N_det_generators
|
||||
if(dress_P(i) <= will_send) sum_f = sum_f + f(i)
|
||||
end do
|
||||
call push_dress_results(zmq_socket_push, -will_send, sum_f, edI_task, edI_index, breve_delta_m, purge_task_id, 1)
|
||||
end if
|
||||
|
||||
!$OMP END SINGLE
|
||||
call disconnect_from_taskserver(zmq_to_qp_run_socket,worker_id)
|
||||
call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket)
|
||||
call end_zmq_push_socket(zmq_socket_push,thread)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do i=0,N_cp+1
|
||||
do i=0,dress_N_cp+1
|
||||
call omp_destroy_lock(lck_sto(i))
|
||||
end do
|
||||
do i=0,comb_teeth+1
|
||||
do i=0,pt2_N_teeth+1
|
||||
call omp_destroy_lock(lck_det(i))
|
||||
end do
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine push_dress_results(zmq_socket_push, ind, cur_cp, delta_loc, int_buf, double_buf, det_buf, N_bufi, task_id)
|
||||
subroutine push_dress_results(zmq_socket_push, m_task, f, edI_task, edI_index, breve_delta_m, task_id, n_tasks)
|
||||
use f77_zmq
|
||||
implicit none
|
||||
|
||||
integer, parameter :: sendt = 4
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_push
|
||||
double precision, intent(inout) :: delta_loc(N_states, N_det, 2)
|
||||
real(kind=4), allocatable :: delta_loc4(:,:,:)
|
||||
double precision, intent(in) :: double_buf(*)
|
||||
integer, intent(in) :: int_buf(*)
|
||||
integer(bit_kind), intent(in) :: det_buf(N_int, 2, *)
|
||||
integer, intent(in) :: N_bufi(3)
|
||||
integer :: N_buf(3)
|
||||
integer, intent(in) :: ind, cur_cp, task_id
|
||||
integer :: rc, i, j, k, l
|
||||
double precision :: contrib(N_states)
|
||||
real(sendt), allocatable :: r4buf(:,:,:)
|
||||
integer, intent(in) :: m_task, f, edI_index(n_tasks)
|
||||
double precision, intent(in) :: breve_delta_m(N_states, N_det, 2), edI_task(n_tasks)
|
||||
integer, intent(in) :: task_id(pt2_n_tasks_max), n_tasks
|
||||
integer :: rc, i, j, k
|
||||
rc = f77_zmq_send( zmq_socket_push, m_task, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push3"
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, ind, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push"
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, cur_cp, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push"
|
||||
|
||||
|
||||
if(cur_cp /= -1) then
|
||||
allocate(r4buf(N_states, N_det, 2))
|
||||
do i=1,2
|
||||
do j=1,N_det
|
||||
do k=1,N_states
|
||||
r4buf(k,j,i) = real(delta_loc(k,j,i), sendt)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, r4buf(1,1,1), sendt*N_states*N_det, ZMQ_SNDMORE)
|
||||
if(rc /= sendt*N_states*N_det) stop "push"
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, r4buf(1,1,2), sendt*N_states*N_det, ZMQ_SNDMORE)
|
||||
if(rc /= sendt*N_states*N_det) stop "push"
|
||||
if(m_task > 0) then
|
||||
rc = f77_zmq_send( zmq_socket_push, n_tasks, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push1"
|
||||
rc = f77_zmq_send( zmq_socket_push, f, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push4"
|
||||
rc = f77_zmq_send( zmq_socket_push, edI_task, 8*n_tasks, ZMQ_SNDMORE)
|
||||
if(rc /= 8*n_tasks) stop "push5"
|
||||
rc = f77_zmq_send( zmq_socket_push, edI_index, 4*n_tasks, 0)
|
||||
if(rc /= 4*n_tasks) stop "push6"
|
||||
else if(m_task == 0) then
|
||||
rc = f77_zmq_send( zmq_socket_push, n_tasks, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push1"
|
||||
rc = f77_zmq_send( zmq_socket_push, task_id, 4*n_tasks, 0)
|
||||
if(rc /= 4*n_tasks) stop "push2"
|
||||
else
|
||||
contrib = 0d0
|
||||
do i=1,N_det
|
||||
contrib(:) += delta_loc(:,i, 1) * psi_coef(i, :)
|
||||
end do
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, contrib, 8*N_states, ZMQ_SNDMORE)
|
||||
if(rc /= 8*N_states) stop "push"
|
||||
|
||||
N_buf = N_bufi
|
||||
!N_buf = (/0,1,0/)
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, N_buf, 4*3, ZMQ_SNDMORE)
|
||||
if(rc /= 4*3) stop "push5"
|
||||
|
||||
if(N_buf(1) > N_dress_int_buffer) stop "run_dress_slave N_buf bad size?"
|
||||
if(N_buf(2) > N_dress_double_buffer) stop "run_dress_slave N_buf bad size?"
|
||||
if(N_buf(3) > N_dress_det_buffer) stop "run_dress_slave N_buf bad size?"
|
||||
|
||||
|
||||
if(N_buf(1) > 0) then
|
||||
rc = f77_zmq_send( zmq_socket_push, int_buf, 4*N_buf(1), ZMQ_SNDMORE)
|
||||
if(rc /= 4*N_buf(1)) stop "push6"
|
||||
end if
|
||||
|
||||
if(N_buf(2) > 0) then
|
||||
rc = f77_zmq_send( zmq_socket_push, double_buf, 8*N_buf(2), ZMQ_SNDMORE)
|
||||
if(rc /= 8*N_buf(2)) stop "push8"
|
||||
end if
|
||||
|
||||
if(N_buf(3) > 0) then
|
||||
rc = f77_zmq_send( zmq_socket_push, det_buf, 2*N_int*bit_kind*N_buf(3), ZMQ_SNDMORE)
|
||||
if(rc /= 2*N_int*bit_kind*N_buf(3)) stop "push10"
|
||||
end if
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_push, f, 4, ZMQ_SNDMORE)
|
||||
if(rc /= 4) stop "push4"
|
||||
rc = f77_zmq_send( zmq_socket_push, breve_delta_m, 8*N_det*N_states*2, ZMQ_SNDMORE)
|
||||
if(rc /= 8*N_det*N_states*2) stop "push6"
|
||||
rc = f77_zmq_send( zmq_socket_push, task_id, 4, 0)
|
||||
if(rc /= 4) stop "push11"
|
||||
if(rc /= 4) stop "push6"
|
||||
end if
|
||||
|
||||
! Activate is zmq_socket_push is a REQ
|
||||
! Activate is zmq_socket_pull is a REP
|
||||
IRP_IF ZMQ_PUSH
|
||||
IRP_ELSE
|
||||
character*(2) :: ok
|
||||
rc = f77_zmq_recv( zmq_socket_push, ok, 2, 0)
|
||||
IRP_ENDIF
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ real(4), real4buf, (N_states, N_det, 2) ]
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine pull_dress_results(zmq_socket_pull, ind, cur_cp, delta_loc, int_buf, double_buf, det_buf, N_buf, task_id, contrib)
|
||||
subroutine pull_dress_results(zmq_socket_pull, m_task, f, edI_task, edI_index, breve_delta_m, task_id, n_tasks)
|
||||
use f77_zmq
|
||||
implicit none
|
||||
integer, parameter :: sendt = 4
|
||||
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
|
||||
integer, intent(out) :: cur_cp
|
||||
double precision, intent(inout) :: delta_loc(N_states, N_det, 2)
|
||||
double precision, intent(out) :: double_buf(*), contrib(N_states)
|
||||
integer, intent(out) :: int_buf(*)
|
||||
integer(bit_kind), intent(out) :: det_buf(N_int, 2, *)
|
||||
integer, intent(out) :: ind
|
||||
integer, intent(out) :: task_id
|
||||
integer, intent(out) :: m_task, f, edI_index(N_det_generators)
|
||||
double precision, intent(out) :: breve_delta_m(N_states, N_det, 2), edI_task(N_det_generators)
|
||||
integer, intent(out) :: task_id(pt2_n_tasks_max), n_tasks
|
||||
integer :: rc, i, j, k
|
||||
integer, intent(out) :: N_buf(3)
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, ind, 4, 0)
|
||||
if(rc /= 4) stop "pulla"
|
||||
rc = f77_zmq_recv( zmq_socket_pull, m_task, 4, 0)
|
||||
if(rc /= 4) stop "pullc"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, cur_cp, 4, 0)
|
||||
if(rc /= 4) stop "pulla"
|
||||
|
||||
|
||||
|
||||
|
||||
if(cur_cp /= -1) then
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, real4buf(1,1,1), N_states*sendt*N_det, 0)
|
||||
if(rc /= sendt*N_states*N_det) stop "pullc"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, real4buf(1,1,2), N_states*sendt*N_det, 0)
|
||||
if(rc /= sendt*N_states*N_det) stop "pulld"
|
||||
|
||||
do i=1,2
|
||||
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,k)
|
||||
do j=1,N_det
|
||||
do k=1,N_states
|
||||
delta_loc(k,j,i) = real(real4buf(k,j,i), 8)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
if(m_task > 0) then
|
||||
rc = f77_zmq_recv( zmq_socket_pull, n_tasks, 4, 0)
|
||||
if(rc /= 4) stop "pullc"
|
||||
rc = f77_zmq_recv( zmq_socket_pull, f, 4, 0)
|
||||
if(rc /= 4) stop "pullc"
|
||||
rc = f77_zmq_recv( zmq_socket_pull, edI_task, 8*n_tasks, 0)
|
||||
if(rc /= 8*n_tasks) stop "pullc"
|
||||
rc = f77_zmq_recv( zmq_socket_pull, edI_index, 4*n_tasks, 0)
|
||||
if(rc /= 4*n_tasks) stop "pullc"
|
||||
else if(m_task==0) then
|
||||
rc = f77_zmq_recv( zmq_socket_pull, n_tasks, 4, 0)
|
||||
if(rc /= 4) stop "pullc"
|
||||
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4*n_tasks, 0)
|
||||
if(rc /= 4*n_tasks) stop "pull4"
|
||||
else
|
||||
rc = f77_zmq_recv( zmq_socket_pull, contrib, 8*N_states, 0)
|
||||
if(rc /= 8*N_states) stop "pullc"
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, N_buf, 4*3, 0)
|
||||
if(rc /= 4*3) stop "pull"
|
||||
if(N_buf(1) > N_dress_int_buffer) stop "run_dress_slave N_buf bad size?"
|
||||
if(N_buf(2) > N_dress_double_buffer) stop "run_dress_slave N_buf bad size?"
|
||||
if(N_buf(3) > N_dress_det_buffer) stop "run_dress_slave N_buf bad size?"
|
||||
|
||||
|
||||
if(N_buf(1) > 0) then
|
||||
rc = f77_zmq_recv( zmq_socket_pull, int_buf, 4*N_buf(1), 0)
|
||||
if(rc /= 4*N_buf(1)) stop "pull1"
|
||||
end if
|
||||
|
||||
if(N_buf(2) > 0) then
|
||||
rc = f77_zmq_recv( zmq_socket_pull, double_buf, 8*N_buf(2), 0)
|
||||
if(rc /= 8*N_buf(2)) stop "pull2"
|
||||
end if
|
||||
|
||||
if(N_buf(3) > 0) then
|
||||
rc = f77_zmq_recv( zmq_socket_pull, det_buf, 2*N_int*bit_kind*N_buf(3), 0)
|
||||
if(rc /= 2*N_int*bit_kind*N_buf(3)) stop "pull3"
|
||||
end if
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull, f, 4, 0)
|
||||
if(rc /= 4) stop "pullc"
|
||||
rc = f77_zmq_recv( zmq_socket_pull, breve_delta_m, 8*N_det*N_states*2, 0)
|
||||
if(rc /= 8*N_det*N_states*2) stop "pullc"
|
||||
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
|
||||
if(rc /= 4) stop "pull4"
|
||||
end if
|
||||
|
||||
@ -26,3 +26,9 @@ doc: Type of zeroth-order Hamiltonian [ EN | Barycentric ]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: EN
|
||||
|
||||
[dress_relative_error]
|
||||
type: Normalized_float
|
||||
doc: Stop stochastic dressing when the relative error is smaller than PT2_relative_error
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.001
|
||||
|
||||
|
||||
@ -1,159 +0,0 @@
|
||||
program shifted_bk
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
double precision, allocatable :: pt2(:)
|
||||
integer :: degree
|
||||
integer :: n_det_before
|
||||
double precision :: threshold_davidson_in
|
||||
|
||||
allocate (pt2(N_states))
|
||||
|
||||
double precision :: hf_energy_ref
|
||||
logical :: has
|
||||
double precision :: relative_error, absolute_error
|
||||
integer :: N_states_p
|
||||
character*(512) :: fmt
|
||||
|
||||
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
|
||||
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
|
||||
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
|
||||
PROVIDE psi_bilinear_matrix_transp_order
|
||||
|
||||
|
||||
pt2 = -huge(1.e0)
|
||||
threshold_davidson_in = threshold_davidson
|
||||
threshold_davidson = threshold_davidson_in * 100.d0
|
||||
SOFT_TOUCH threshold_davidson
|
||||
|
||||
call diagonalize_CI_dressed
|
||||
call save_wavefunction
|
||||
|
||||
call ezfio_has_hartree_fock_energy(has)
|
||||
if (has) then
|
||||
call ezfio_get_hartree_fock_energy(hf_energy_ref)
|
||||
else
|
||||
hf_energy_ref = ref_bitmask_energy
|
||||
endif
|
||||
|
||||
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_dressed
|
||||
call save_wavefunction
|
||||
N_states_p = min(N_det,N_states)
|
||||
endif
|
||||
|
||||
n_det_before = 0
|
||||
|
||||
character*(8) :: pt2_string
|
||||
double precision :: threshold_selectors_save, threshold_generators_save
|
||||
threshold_selectors_save = threshold_selectors
|
||||
threshold_generators_save = threshold_generators
|
||||
double precision :: error(N_states), energy(N_states)
|
||||
error = 0.d0
|
||||
|
||||
threshold_selectors = 1.d0
|
||||
threshold_generators = 1d0
|
||||
|
||||
if (.True.) then
|
||||
pt2_string = '(sh-Bk) '
|
||||
do while ( (N_det < N_det_max) )
|
||||
write(*,'(A)') '--------------------------------------------------------------------------------'
|
||||
|
||||
N_det_delta_ij = N_det
|
||||
|
||||
do i=1,N_states
|
||||
energy(i) = psi_energy(i)+nuclear_repulsion
|
||||
enddo
|
||||
|
||||
PROVIDE delta_ij_tmp
|
||||
call delta_ij_done()
|
||||
|
||||
call diagonalize_ci_dressed
|
||||
do i=1,N_states
|
||||
pt2(i) = ci_energy_dressed(i) - energy(i)
|
||||
enddo
|
||||
|
||||
N_states_p = min(N_det,N_states)
|
||||
|
||||
print *, ''
|
||||
print '(A,I12)', 'Summary at N_det = ', N_det
|
||||
print '(A)', '-----------------------------------'
|
||||
print *, ''
|
||||
print *, ''
|
||||
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))'
|
||||
write(*,fmt) ('State',k, k=1,N_states_p)
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))'
|
||||
write(*,fmt) '# E ', energy(1:N_states_p)
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', energy(1:N_states_p)-energy(1)
|
||||
write(*,fmt) '# Excit. (eV)', (energy(1:N_states_p)-energy(1))*27.211396641308d0
|
||||
endif
|
||||
write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
|
||||
write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
|
||||
write(*,'(A)') '#'
|
||||
write(*,fmt) '# E+PT2 ', (energy(k)+pt2(k),error(k), k=1,N_states_p)
|
||||
if (N_states_p > 1) then
|
||||
write(*,fmt) '# Excit. (au)', ( (energy(k)+pt2(k)-energy(1)-pt2(1)), &
|
||||
dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
|
||||
write(*,fmt) '# Excit. (eV)', ( (energy(k)+pt2(k)-energy(1)-pt2(1))*27.211396641308d0, &
|
||||
dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
|
||||
endif
|
||||
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
|
||||
write(*,fmt)
|
||||
print *, ''
|
||||
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
|
||||
do k=1, N_states_p
|
||||
print*,'State ',k
|
||||
print *, 'PT2 = ', pt2(k)
|
||||
print *, 'E = ', energy(k)
|
||||
print *, 'E+PT2'//pt2_string//' = ', energy(k)+pt2(k)
|
||||
enddo
|
||||
|
||||
print *, '-----'
|
||||
if(N_states.gt.1)then
|
||||
print *, 'Variational Energy difference (au | eV)'
|
||||
do i=2, N_states_p
|
||||
print*,'Delta E = ', (energy(i) - energy(1)), &
|
||||
(energy(i) - energy(1)) * 27.211396641308d0
|
||||
enddo
|
||||
print *, '-----'
|
||||
print*, 'Variational + perturbative Energy difference (au | eV)'
|
||||
do i=2, N_states_p
|
||||
print*,'Delta E = ', (energy(i)+ pt2(i) - (energy(1) + pt2(1))), &
|
||||
(energy(i)+ pt2(i) - (energy(1) + pt2(1))) * 27.211396641308d0
|
||||
enddo
|
||||
endif
|
||||
call ezfio_set_shiftedbk_energy_pt2(energy(1)+pt2(1))
|
||||
! call dump_fci_iterations_value(N_det,energy,pt2)
|
||||
|
||||
n_det_before = N_det
|
||||
|
||||
PROVIDE psi_coef
|
||||
PROVIDE psi_det
|
||||
PROVIDE psi_det_sorted
|
||||
|
||||
if (N_det >= N_det_max) then
|
||||
threshold_davidson = threshold_davidson_in
|
||||
end if
|
||||
call save_wavefunction
|
||||
call ezfio_set_shiftedbk_energy(energy(1))
|
||||
call ezfio_set_shiftedbk_energy_pt2(ci_energy_dressed(1))
|
||||
enddo
|
||||
endif
|
||||
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
@ -300,7 +300,6 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili
|
||||
haa = diag_H_mat_elem_fock(psi_det_generators(1,1,i_gen),alpha,fock_diag_tmp_(1,1,iproc),N_int)
|
||||
|
||||
call dress_with_alpha_(Nstates, Ndet, Nint, delta_ij_loc, minilist, det_minilist, n_minilist, alpha, haa, contrib, c_alpha, iproc)
|
||||
|
||||
slave_sum_alpha2(:,iproc) += c_alpha(:)**2
|
||||
if(contrib < sb(iproc)%mini) then
|
||||
call add_to_selection_buffer(sb(iproc), alpha, contrib)
|
||||
|
||||
@ -8,7 +8,7 @@ default: 1.e-12
|
||||
type: States_number
|
||||
doc: Number of states to consider during the Davdison diagonalization
|
||||
default: 4
|
||||
interface: ezfio,provider,ocaml
|
||||
interface: ezfio,ocaml
|
||||
|
||||
[davidson_sze_max]
|
||||
type: Strictly_positive_int
|
||||
|
||||
35
src/Davidson/ezfio.irp.f
Normal file
35
src/Davidson/ezfio.irp.f
Normal file
@ -0,0 +1,35 @@
|
||||
BEGIN_PROVIDER [ integer, n_states_diag ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of states to consider during the Davdison diagonalization
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
if (mpi_master) then
|
||||
|
||||
call ezfio_has_davidson_n_states_diag(has)
|
||||
if (has) then
|
||||
call ezfio_get_davidson_n_states_diag(n_states_diag)
|
||||
else
|
||||
print *, 'davidson/n_states_diag not found in EZFIO file'
|
||||
stop 1
|
||||
endif
|
||||
n_states_diag = max(N_states, N_states_diag)
|
||||
endif
|
||||
IRP_IF MPI
|
||||
include 'mpif.h'
|
||||
integer :: ierr
|
||||
call MPI_BCAST( n_states_diag, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
|
||||
if (ierr /= MPI_SUCCESS) then
|
||||
stop 'Unable to read n_states_diag with MPI'
|
||||
endif
|
||||
IRP_ENDIF
|
||||
|
||||
call write_time(6)
|
||||
if (mpi_master) then
|
||||
write(6, *) 'Read n_states_diag'
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -28,6 +28,7 @@ subroutine H_S2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze)
|
||||
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
|
||||
allocate(u_t(N_st,N_det),v_t(N_st,N_det),s_t(N_st,N_det))
|
||||
|
||||
do k=1,N_st
|
||||
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
|
||||
enddo
|
||||
|
||||
@ -167,8 +167,7 @@ subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
end select
|
||||
end
|
||||
|
||||
|
||||
subroutine get_double_excitation(det1,det2,exc,phase,Nint)
|
||||
subroutine get_double_excitation_ref(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -312,6 +311,137 @@ subroutine get_double_excitation(det1,det2,exc,phase,Nint)
|
||||
|
||||
end
|
||||
|
||||
subroutine get_phasemask_bit(det1, pm, Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint,2)
|
||||
integer(bit_kind), intent(out) :: pm(Nint,2)
|
||||
integer(bit_kind) :: tmp
|
||||
integer :: ispin, i
|
||||
do ispin=1,2
|
||||
tmp = 0_8
|
||||
do i=1,Nint
|
||||
pm(i,ispin) = xor(det1(i,ispin), ishft(det1(i,ispin), 1))
|
||||
pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 2))
|
||||
pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 4))
|
||||
pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 8))
|
||||
pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 16))
|
||||
pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 32))
|
||||
pm(i,ispin) = xor(pm(i,ispin), tmp)
|
||||
if(iand(popcnt(det1(i,ispin)), 1) == 1) tmp = not(tmp)
|
||||
end do
|
||||
end do
|
||||
end subroutine
|
||||
|
||||
|
||||
subroutine get_double_excitation(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Returns the two excitation operators between two doubly excited determinants and the phase
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: det1(Nint,2)
|
||||
integer(bit_kind), intent(in) :: det2(Nint,2)
|
||||
integer, intent(out) :: exc(0:2,2,2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: tz
|
||||
integer :: l, ispin, idx_hole, idx_particle, ishift
|
||||
integer :: nperm
|
||||
integer :: i,j,k,m,n
|
||||
integer :: high, low
|
||||
integer :: a,b,c,d
|
||||
integer(bit_kind) :: hole, particle, tmp, pm(Nint,2)
|
||||
double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /)
|
||||
double precision :: refaz
|
||||
logical :: ok
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
|
||||
!do ispin=1,2
|
||||
!tmp = 0_8
|
||||
!do i=1,Nint
|
||||
! pm(i,ispin) = xor(det1(i,ispin), ishft(det1(i,ispin), 1))
|
||||
! pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 2))
|
||||
! pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 4))
|
||||
! pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 8))
|
||||
! pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 16))
|
||||
! pm(i,ispin) = xor(pm(i,ispin), ishft(pm(i,ispin), 32))
|
||||
! pm(i,ispin) = xor(pm(i,ispin), tmp)
|
||||
! if(iand(popcnt(det1(i,ispin)), 1) == 1) tmp = not(tmp)
|
||||
!end do
|
||||
!end do
|
||||
call get_phasemask_bit(det1, pm, Nint)
|
||||
nperm = 0
|
||||
exc(0,1,1) = 0
|
||||
exc(0,2,1) = 0
|
||||
exc(0,1,2) = 0
|
||||
exc(0,2,2) = 0
|
||||
do ispin = 1,2
|
||||
idx_particle = 0
|
||||
idx_hole = 0
|
||||
ishift = 1-bit_kind_size
|
||||
!par = 0
|
||||
do l=1,Nint
|
||||
ishift = ishift + bit_kind_size
|
||||
if (det1(l,ispin) == det2(l,ispin)) then
|
||||
cycle
|
||||
endif
|
||||
tmp = xor( det1(l,ispin), det2(l,ispin) )
|
||||
particle = iand(tmp, det2(l,ispin))
|
||||
hole = iand(tmp, det1(l,ispin))
|
||||
do while (particle /= 0_bit_kind)
|
||||
tz = trailz(particle)
|
||||
nperm = nperm + iand(ishft(pm(l,ispin), -tz), 1)
|
||||
idx_particle = idx_particle + 1
|
||||
exc(0,2,ispin) = exc(0,2,ispin) + 1
|
||||
exc(idx_particle,2,ispin) = tz+ishift
|
||||
particle = iand(particle,particle-1_bit_kind)
|
||||
enddo
|
||||
if (iand(exc(0,1,ispin),exc(0,2,ispin))==2) then ! exc(0,1,ispin)==2 or exc(0,2,ispin)==2
|
||||
exit
|
||||
endif
|
||||
do while (hole /= 0_bit_kind)
|
||||
tz = trailz(hole)
|
||||
nperm = nperm + iand(ishft(pm(l,ispin), -tz), 1)
|
||||
idx_hole = idx_hole + 1
|
||||
exc(0,1,ispin) = exc(0,1,ispin) + 1
|
||||
exc(idx_hole,1,ispin) = tz+ishift
|
||||
hole = iand(hole,hole-1_bit_kind)
|
||||
enddo
|
||||
if (iand(exc(0,1,ispin),exc(0,2,ispin))==2) then ! exc(0,1,ispin)==2 or exc(0,2,ispin)
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
select case (exc(0,1,ispin))
|
||||
case(0)
|
||||
cycle
|
||||
|
||||
case(1)
|
||||
if(exc(1,1,ispin) < exc(1,2,ispin)) nperm = nperm+1
|
||||
|
||||
case (2)
|
||||
a = exc(1,1,ispin)
|
||||
b = exc(1,2,ispin)
|
||||
c = exc(2,1,ispin)
|
||||
d = exc(2,2,ispin)
|
||||
|
||||
if(min(a,c) > max(b,d) .or. min(b,d) > max(a,c) .or. (a-b)*(c-d)<0) then
|
||||
nperm = nperm + 1
|
||||
end if
|
||||
exit
|
||||
end select
|
||||
|
||||
enddo
|
||||
phase = phase_dble(iand(nperm,1))
|
||||
!call get_double_excitation_ref(det1,det2,exc,refaz,Nint)
|
||||
!if(phase == refaz) then
|
||||
! print *, "phase", phase, refaz, n, exc(0,1,1)
|
||||
!end if
|
||||
end
|
||||
|
||||
subroutine get_mono_excitation(det1,det2,exc,phase,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
Loading…
Reference in New Issue
Block a user