mirror of
https://github.com/LCPQ/quantum_package
synced 2024-08-24 21:41:46 +02:00
Save wavefunction in CISD
This commit is contained in:
commit
199d7ae0d8
@ -15,6 +15,7 @@ function _install()
|
||||
make -j 8 || exit 1
|
||||
mv libf77zmq.a "${QP_ROOT}"/lib || exit 1
|
||||
mv libf77zmq.so "${QP_ROOT}"/lib || exit 1
|
||||
cp f77_zmq.h "${QP_ROOT}"/src/ZMQ/
|
||||
cd -
|
||||
return 0
|
||||
}
|
||||
|
10
plugins/CASSCF/EZFIO.cfg
Normal file
10
plugins/CASSCF/EZFIO.cfg
Normal file
@ -0,0 +1,10 @@
|
||||
[energy]
|
||||
type: double precision
|
||||
doc: "Calculated CAS-SCF energy"
|
||||
interface: ezfio
|
||||
|
||||
[energy_pt2]
|
||||
type: double precision
|
||||
doc: "Calculated selected CAS-SCF energy with PT2 correction"
|
||||
interface: ezfio
|
||||
|
39
plugins/CASSCF/H_apply.irp.f
Normal file
39
plugins/CASSCF/H_apply.irp.f
Normal file
@ -0,0 +1,39 @@
|
||||
use bitmasks
|
||||
BEGIN_SHELL [ /usr/bin/env python ]
|
||||
from generate_h_apply import *
|
||||
|
||||
s = H_apply("CAS_SD")
|
||||
print s
|
||||
|
||||
s = H_apply("CAS_SD_selected_no_skip")
|
||||
s.set_selection_pt2("epstein_nesbet_2x2")
|
||||
s.unset_skip()
|
||||
print s
|
||||
|
||||
s = H_apply("CAS_SD_selected")
|
||||
s.set_selection_pt2("epstein_nesbet_2x2")
|
||||
print s
|
||||
|
||||
s = H_apply("CAS_SD_PT2")
|
||||
s.set_perturbation("epstein_nesbet_2x2")
|
||||
print s
|
||||
|
||||
|
||||
s = H_apply("CAS_S",do_double_exc=False)
|
||||
print s
|
||||
|
||||
s = H_apply("CAS_S_selected_no_skip",do_double_exc=False)
|
||||
s.set_selection_pt2("epstein_nesbet_2x2")
|
||||
s.unset_skip()
|
||||
print s
|
||||
|
||||
s = H_apply("CAS_S_selected",do_double_exc=False)
|
||||
s.set_selection_pt2("epstein_nesbet_2x2")
|
||||
print s
|
||||
|
||||
s = H_apply("CAS_S_PT2",do_double_exc=False)
|
||||
s.set_perturbation("epstein_nesbet_2x2")
|
||||
print s
|
||||
|
||||
END_SHELL
|
||||
|
1
plugins/CASSCF/NEEDED_CHILDREN_MODULES
Normal file
1
plugins/CASSCF/NEEDED_CHILDREN_MODULES
Normal file
@ -0,0 +1 @@
|
||||
Generators_CAS Perturbation Selectors_full
|
20
plugins/CASSCF/README.rst
Normal file
20
plugins/CASSCF/README.rst
Normal file
@ -0,0 +1,20 @@
|
||||
======
|
||||
CASSCF
|
||||
======
|
||||
|
||||
This module is not a "real" CAS-SCF. It is an orbital optimization step done by :
|
||||
|
||||
1) Doing the CAS+SD
|
||||
2) Taking one-electron density matrix
|
||||
3) Cancelling all active-active rotations
|
||||
4) Finding the order which matches with the input MOs
|
||||
|
||||
|
||||
Needed Modules
|
||||
==============
|
||||
.. Do not edit this section It was auto-generated
|
||||
.. by the `update_README.py` script.
|
||||
Documentation
|
||||
=============
|
||||
.. Do not edit this section It was auto-generated
|
||||
.. by the `update_README.py` script.
|
220
plugins/CASSCF/casscf.irp.f
Normal file
220
plugins/CASSCF/casscf.irp.f
Normal file
@ -0,0 +1,220 @@
|
||||
program casscf
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Optimize MOs and CI coefficients of the CAS
|
||||
END_DOC
|
||||
|
||||
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
|
||||
integer(bit_kind), allocatable :: generators_bitmask_save(:,:,:,:)
|
||||
|
||||
integer :: degree, N_generators_bitmask_save, N_det_ci
|
||||
double precision :: E_old, E_CI
|
||||
double precision :: selection_criterion_save, selection_criterion_min_save
|
||||
|
||||
integer :: N_det_old
|
||||
integer :: i, j, k, l
|
||||
integer :: i_bit, j_bit, i_int, j_int
|
||||
integer(bit_kind), allocatable :: bit_tmp(:), cas_bm(:)
|
||||
character*(64) :: label
|
||||
|
||||
allocate( pt2(N_states), norm_pert(N_states),H_pert_diag(N_states) )
|
||||
allocate( generators_bitmask_save(N_int,2,6,N_generators_bitmask) )
|
||||
allocate( bit_tmp(N_int), cas_bm(N_int) )
|
||||
|
||||
PROVIDE N_det_cas
|
||||
N_det_old = 0
|
||||
pt2 = 1.d0
|
||||
E_CI = 1.d0
|
||||
E_old = 0.d0
|
||||
diag_algorithm = "Lapack"
|
||||
selection_criterion_save = selection_criterion
|
||||
selection_criterion_min_save = selection_criterion_min
|
||||
|
||||
|
||||
cas_bm = 0_bit_kind
|
||||
do i=1,N_cas_bitmask
|
||||
do j=1,N_int
|
||||
cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,1,i))
|
||||
cas_bm(j) = ior(cas_bm(j), cas_bitmask(j,2,i))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Save CAS-SD bitmask
|
||||
generators_bitmask_save = generators_bitmask
|
||||
N_generators_bitmask_save = N_generators_bitmask
|
||||
|
||||
! Set the CAS bitmask
|
||||
do i=1,6
|
||||
generators_bitmask(:,:,i,:) = cas_bitmask
|
||||
enddo
|
||||
N_generators_bitmask = N_cas_bitmask
|
||||
SOFT_TOUCH generators_bitmask N_generators_bitmask
|
||||
|
||||
|
||||
! If the number of dets already in the file is larger than the requested
|
||||
! number of determinants, truncate the wf
|
||||
if (N_det > N_det_max) then
|
||||
call diagonalize_CI
|
||||
call save_wavefunction
|
||||
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
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
print *, 'PT2 = ', pt2
|
||||
print *, 'E = ', CI_energy
|
||||
print *, 'E+PT2 = ', CI_energy+pt2
|
||||
print *, '-----'
|
||||
endif
|
||||
|
||||
! Start MCSCF iteration
|
||||
|
||||
! CAS-CI
|
||||
! ------
|
||||
|
||||
E_old = E_CI
|
||||
|
||||
! Reset the selection criterion
|
||||
selection_criterion = selection_criterion_save
|
||||
selection_criterion_min = selection_criterion_min_save
|
||||
SOFT_TOUCH selection_criterion_min selection_criterion selection_criterion_factor
|
||||
|
||||
! Set the CAS bitmask
|
||||
do i=1,6
|
||||
generators_bitmask(:,:,i,:) = cas_bitmask
|
||||
enddo
|
||||
N_generators_bitmask = N_cas_bitmask
|
||||
SOFT_TOUCH generators_bitmask N_generators_bitmask
|
||||
|
||||
do while (N_det < N_det_max.and.maxval(abs(pt2(1:N_states))) > pt2_max)
|
||||
N_det_old = N_det
|
||||
call H_apply_CAS_SD_selected_no_skip(pt2, norm_pert, H_pert_diag, N_states)
|
||||
|
||||
PROVIDE psi_coef
|
||||
PROVIDE psi_det
|
||||
PROVIDE psi_det_sorted
|
||||
|
||||
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
|
||||
endif
|
||||
call diagonalize_CI
|
||||
call save_wavefunction
|
||||
print *, '======'
|
||||
print *, 'CAS-CI'
|
||||
print *, '======'
|
||||
print *, ''
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
print *, 'PT2 = ', pt2
|
||||
print *, 'E(CAS) = ', CI_energy
|
||||
print *, 'E(CAS)+PT2 = ', CI_energy+pt2
|
||||
print *, '-----'
|
||||
print *, ''
|
||||
E_CI = sum(CI_energy(1:N_states)+pt2(1:N_states))/dble(N_states)
|
||||
|
||||
call ezfio_set_casscf_energy(CI_energy(1))
|
||||
if (abort_all) then
|
||||
exit
|
||||
endif
|
||||
if (N_det == N_det_old) then
|
||||
exit
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
! Super-CI
|
||||
! --------
|
||||
|
||||
selection_criterion_min = 1.d-12
|
||||
selection_criterion = 1.d-12
|
||||
|
||||
! Set the CAS bitmask
|
||||
generators_bitmask = generators_bitmask_save
|
||||
N_generators_bitmask = N_generators_bitmask_save
|
||||
SOFT_TOUCH generators_bitmask N_generators_bitmask selection_criterion selection_criterion_min selection_criterion_factor
|
||||
|
||||
N_det_ci = N_det
|
||||
|
||||
call H_apply_CAS_SD_selected(pt2, norm_pert, H_pert_diag, N_states)
|
||||
|
||||
|
||||
do i=1,mo_tot_num
|
||||
i_int = ishft(i-1,-bit_kind_shift)+1
|
||||
i_bit = j-ishft(i_int-1,bit_kind_shift)-1
|
||||
bit_tmp(:) = 0_bit_kind
|
||||
bit_tmp(i_int) = ibset(0_bit_kind,i_bit)
|
||||
if (iand(bit_tmp(i_int), cas_bm(i_int)) == 0_bit_kind) then
|
||||
! Not a CAS MO
|
||||
cycle
|
||||
endif
|
||||
|
||||
do j=1,mo_tot_num
|
||||
if (j == i) then
|
||||
cycle
|
||||
endif
|
||||
j_int = ishft(j-1,-bit_kind_shift)+1
|
||||
j_bit = j-ishft(j_int-1,bit_kind_shift)-1
|
||||
bit_tmp(:) = 0_bit_kind
|
||||
bit_tmp(j_int) = ibset(0_bit_kind,j_bit)
|
||||
if (iand(bit_tmp(j_int), cas_bm(j_int)) == 0_bit_kind) then
|
||||
! Not a CAS MO
|
||||
cycle
|
||||
endif
|
||||
! Now, both i and j are MOs of the CAS. De-couple them in the DM
|
||||
one_body_dm_mo(i,j) = 0.d0
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
SOFT_TOUCH one_body_dm_mo
|
||||
|
||||
double precision :: mx, ov
|
||||
double precision, allocatable :: mo_coef_old(:,:)
|
||||
integer, allocatable :: iorder(:)
|
||||
logical, allocatable :: selected(:)
|
||||
allocate( mo_coef_old(size(mo_coef,1), size(mo_coef,2)), iorder(mo_tot_num), selected(mo_tot_num) )
|
||||
mo_coef_old = mo_coef
|
||||
label = "Canonical"
|
||||
call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,-1)
|
||||
selected = .False.
|
||||
do j=1,mo_tot_num
|
||||
mx = -1.d0
|
||||
iorder(j) = j
|
||||
do i=1,mo_tot_num
|
||||
if (selected(i)) then
|
||||
cycle
|
||||
endif
|
||||
ov = 0.d0
|
||||
do l=1,ao_num
|
||||
do k=1,ao_num
|
||||
ov = ov + mo_coef_old(k,j) * ao_overlap(k,l) * mo_coef(l,i)
|
||||
enddo
|
||||
enddo
|
||||
ov= dabs(ov)
|
||||
if (ov > mx) then
|
||||
mx = ov
|
||||
iorder(j) = i
|
||||
endif
|
||||
enddo
|
||||
selected( iorder(j) ) = .True.
|
||||
enddo
|
||||
mo_coef_old = mo_coef
|
||||
do i=1,mo_tot_num
|
||||
mo_coef(:,i) = mo_coef_old(:,iorder(i))
|
||||
enddo
|
||||
|
||||
call save_mos
|
||||
|
||||
call write_double(6,E_CI,"Energy(CAS)")
|
||||
|
||||
deallocate( mo_coef_old )
|
||||
deallocate( pt2, norm_pert,H_pert_diag )
|
||||
deallocate( generators_bitmask_save )
|
||||
deallocate( bit_tmp, cas_bm, iorder )
|
||||
end
|
@ -41,8 +41,8 @@ program cisd
|
||||
N_det = min(N_det,N_det_max)
|
||||
touch N_det psi_det psi_coef
|
||||
call diagonalize_CI
|
||||
deallocate(pt2,norm_pert,H_pert_diag)
|
||||
call save_wavefunction
|
||||
call ezfio_set_cisd_selected_energy(CI_energy)
|
||||
call ezfio_set_cisd_selected_energy_pt2(CI_energy+pt2)
|
||||
deallocate(pt2,norm_pert,H_pert_diag)
|
||||
end
|
||||
|
@ -1,6 +1,7 @@
|
||||
program save_for_qmc
|
||||
read_wf = .True.
|
||||
TOUCH read_wf
|
||||
print *, "N_det = ", N_det
|
||||
call write_spindeterminants
|
||||
if (do_pseudo) then
|
||||
call write_pseudopotential
|
||||
|
@ -36,8 +36,8 @@ except ImportError:
|
||||
from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
|
||||
|
||||
LIB = "" # join(QP_ROOT, "lib", "rdtsc.o")
|
||||
EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a")
|
||||
ZMQ_LIB = join(QP_ROOT, "lib", "libzmq.a") + " " + join(QP_ROOT, "lib", "libf77zmq.a")
|
||||
EZFIO_LIB = join(QP_ROOT, "lib", "libezfio_irp.a")
|
||||
ZMQ_LIB = join(QP_ROOT, "lib", "libf77zmq.a") + " " + join(QP_ROOT, "lib", "libzmq.a") + " -lstdc++ -lrt"
|
||||
ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja")
|
||||
|
||||
header = r"""#
|
||||
@ -262,7 +262,7 @@ def ninja_ezfio_rule():
|
||||
l_flag = ["export {0}='${0}'".format(flag)
|
||||
for flag in ["FC", "FCFLAGS", "IRPF90"]]
|
||||
|
||||
install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio.a")
|
||||
install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio_irp.a")
|
||||
l_cmd = ["cd {0}".format(QP_EZFIO)] + l_flag
|
||||
l_cmd += ["rm -f make.config ; ninja && ln -sf {0} {1}".format(install_lib_ezfio, EZFIO_LIB)]
|
||||
|
||||
|
@ -88,7 +88,7 @@ def get_l_module_descendant(d_child, l_module):
|
||||
except KeyError:
|
||||
print >> sys.stderr, "Error: "
|
||||
print >> sys.stderr, "`{0}` is not a submodule".format(module)
|
||||
print >> sys.stderr, "Check the typo (orthograph, case, '/', etc.) "
|
||||
print >> sys.stderr, "Check the typo (spelling, case, '/', etc.) "
|
||||
sys.exit(1)
|
||||
|
||||
return list(set(l))
|
||||
|
@ -262,13 +262,7 @@ END_PROVIDER
|
||||
logical :: exists
|
||||
integer :: j,i
|
||||
integer :: i_hole,i_part,i_gen
|
||||
PROVIDE ezfio_filename
|
||||
!do j = 1, N_int
|
||||
! inact_bitmask(j,1) = xor(generators_bitmask(j,1,1,1),cas_bitmask(j,1,1))
|
||||
! inact_bitmask(j,2) = xor(generators_bitmask(j,2,1,1),cas_bitmask(j,2,1))
|
||||
! virt_bitmask(j,1) = xor(generators_bitmask(j,1,2,1),cas_bitmask(j,1,1))
|
||||
! virt_bitmask(j,2) = xor(generators_bitmask(j,2,2,1),cas_bitmask(j,2,1))
|
||||
!enddo
|
||||
|
||||
n_inact_orb = 0
|
||||
n_virt_orb = 0
|
||||
if(N_generators_bitmask == 1)then
|
||||
|
@ -1 +1 @@
|
||||
Pseudo Bitmask
|
||||
Pseudo Bitmask ZMQ
|
||||
|
@ -301,7 +301,7 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value)
|
||||
double precision :: thresh
|
||||
thresh = ao_integrals_threshold
|
||||
|
||||
integer :: n_centers, i
|
||||
integer :: i
|
||||
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
buffer_value = 0._integral_kind
|
||||
@ -329,6 +329,7 @@ end
|
||||
|
||||
BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
|
||||
implicit none
|
||||
use f77_zmq
|
||||
use map_module
|
||||
BEGIN_DOC
|
||||
! Map of Atomic integrals
|
||||
@ -345,9 +346,8 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
|
||||
integer(key_kind),allocatable :: buffer_i(:)
|
||||
integer,parameter :: size_buffer = 1024*64
|
||||
real(integral_kind),allocatable :: buffer_value(:)
|
||||
integer(omp_lock_kind) :: lock
|
||||
|
||||
integer :: n_integrals, n_centers, thread_num
|
||||
integer :: n_integrals, rc
|
||||
integer :: jl_pairs(2,ao_num*(ao_num+1)/2), kk, m, j1, i1, lmax
|
||||
|
||||
integral = ao_bielec_integral(1,1,1,1)
|
||||
@ -363,120 +363,61 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
|
||||
endif
|
||||
endif
|
||||
|
||||
kk=1
|
||||
do l = 1, ao_num ! r2
|
||||
do j = 1, l ! r2
|
||||
jl_pairs(1,kk) = j
|
||||
jl_pairs(2,kk) = l
|
||||
kk += 1
|
||||
enddo
|
||||
enddo
|
||||
|
||||
PROVIDE progress_bar
|
||||
call omp_init_lock(lock)
|
||||
lmax = ao_num*(ao_num+1)/2
|
||||
print*, 'Providing the AO integrals'
|
||||
call wall_time(wall_0)
|
||||
call wall_time(wall_1)
|
||||
call cpu_time(cpu_1)
|
||||
call start_progress(lmax,'AO integrals (MB)',0.d0)
|
||||
!$OMP PARALLEL PRIVATE(i,j,k,l,kk, &
|
||||
!$OMP integral,buffer_i,buffer_value,n_integrals, &
|
||||
!$OMP cpu_2,wall_2,i1,j1,thread_num) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, &
|
||||
!$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, &
|
||||
!$OMP ao_overlap_abs,ao_overlap,abort_here, &
|
||||
!$OMP wall_0,progress_bar,progress_value, &
|
||||
!$OMP ao_bielec_integral_schwartz)
|
||||
|
||||
allocate(buffer_i(size_buffer))
|
||||
allocate(buffer_value(size_buffer))
|
||||
n_integrals = 0
|
||||
!$ thread_num = omp_get_thread_num()
|
||||
integer(ZMQ_PTR) :: zmq_socket_rep_inproc, zmq_socket_push_inproc
|
||||
zmq_socket_rep_inproc = f77_zmq_socket(zmq_context, ZMQ_REP)
|
||||
rc = f77_zmq_bind(zmq_socket_rep_inproc, 'inproc://req_rep')
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_socket_rep_inproc'
|
||||
endif
|
||||
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do kk=1,lmax
|
||||
IRP_IF COARRAY
|
||||
if (mod(kk-this_image(),num_images()) /= 0) then
|
||||
cycle
|
||||
endif
|
||||
IRP_ENDIF
|
||||
if (abort_here) then
|
||||
cycle
|
||||
endif
|
||||
if (thread_num == 0) then
|
||||
progress_bar(1) = kk
|
||||
endif
|
||||
j = jl_pairs(1,kk)
|
||||
l = jl_pairs(2,kk)
|
||||
j1 = j+ishft(l*l-l,-1)
|
||||
if (ao_overlap_abs(j,l) < thresh) then
|
||||
cycle
|
||||
endif
|
||||
do k = 1, ao_num ! r1
|
||||
i1 = ishft(k*k-k,-1)
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
do i = 1, k
|
||||
i1 += 1
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
|
||||
cycle
|
||||
endif
|
||||
if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then
|
||||
cycle
|
||||
endif
|
||||
!DIR$ FORCEINLINE
|
||||
integral = ao_bielec_integral(i,k,j,l)
|
||||
if (abs(integral) < thresh) then
|
||||
cycle
|
||||
endif
|
||||
n_integrals += 1
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
|
||||
buffer_value(n_integrals) = integral
|
||||
if (n_integrals > 1024 ) then
|
||||
if (omp_test_lock(lock)) then
|
||||
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
|
||||
n_integrals = 0
|
||||
call omp_unset_lock(lock)
|
||||
endif
|
||||
endif
|
||||
if (n_integrals == size(buffer_i)) then
|
||||
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
|
||||
n_integrals = 0
|
||||
endif
|
||||
integer(ZMQ_PTR) :: thread(0:nproc)
|
||||
external :: ao_bielec_integrals_in_map_slave, ao_bielec_integrals_in_map_collector
|
||||
rc = pthread_create( thread(0), ao_bielec_integrals_in_map_collector )
|
||||
! Create client threads
|
||||
do i=1,nproc
|
||||
rc = pthread_create( thread(i), ao_bielec_integrals_in_map_slave )
|
||||
enddo
|
||||
enddo
|
||||
call wall_time(wall_2)
|
||||
|
||||
if (thread_num == 0) then
|
||||
if (wall_2 - wall_0 > 1.d0) then
|
||||
wall_0 = wall_2
|
||||
print*, 100.*float(kk)/float(lmax), '% in ', &
|
||||
wall_2-wall_1, 's', map_mb(ao_integrals_map) ,'MB'
|
||||
progress_value = dble(map_mb(ao_integrals_map))
|
||||
endif
|
||||
endif
|
||||
character*(64) :: message_string
|
||||
|
||||
do l = ao_num, 1, -1
|
||||
rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0)
|
||||
print *, l
|
||||
! TODO : error handling
|
||||
ASSERT (rc >= 0)
|
||||
ASSERT (message == 'get_ao_integrals')
|
||||
rc = f77_zmq_send( zmq_socket_rep_inproc, l, 4, 0)
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
|
||||
deallocate(buffer_i)
|
||||
deallocate(buffer_value)
|
||||
!$OMP END PARALLEL
|
||||
call omp_destroy_lock(lock)
|
||||
call stop_progress
|
||||
if (abort_here) then
|
||||
stop 'Aborting in AO integrals calculation'
|
||||
do i=1,nproc
|
||||
rc = f77_zmq_recv( zmq_socket_rep_inproc, message_string, 64, 0)
|
||||
! TODO : error handling
|
||||
ASSERT (rc >= 0)
|
||||
ASSERT (message == 'get_ao_integrals')
|
||||
rc = f77_zmq_send( zmq_socket_rep_inproc, 0, 4, 0)
|
||||
enddo
|
||||
! TODO terminer thread(0)
|
||||
|
||||
rc = f77_zmq_unbind(zmq_socket_rep_inproc, 'inproc://req_rep')
|
||||
do i=1,nproc
|
||||
rc = pthread_join( thread(i) )
|
||||
enddo
|
||||
|
||||
zmq_socket_push_inproc = f77_zmq_socket(zmq_context, ZMQ_PUSH)
|
||||
rc = f77_zmq_connect(zmq_socket_push_inproc, 'inproc://push_pull')
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_socket_push_inproc'
|
||||
endif
|
||||
IRP_IF COARRAY
|
||||
print*, 'Communicating the map'
|
||||
call communicate_ao_integrals()
|
||||
IRP_ENDIF COARRAY
|
||||
rc = f77_zmq_send( zmq_socket_push_inproc, -1, 4, ZMQ_SNDMORE)
|
||||
rc = f77_zmq_send( zmq_socket_push_inproc, 0_key_kind, key_kind, ZMQ_SNDMORE)
|
||||
rc = f77_zmq_send( zmq_socket_push_inproc, 0_integral_kind, integral_kind, 0)
|
||||
|
||||
rc = pthread_join( thread(0) )
|
||||
|
||||
print*, 'Sorting the map'
|
||||
call map_sort(ao_integrals_map)
|
||||
call cpu_time(cpu_2)
|
||||
@ -1256,3 +1197,57 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
|
||||
implicit none
|
||||
use map_module
|
||||
BEGIN_DOC
|
||||
! Parallel client for AO integrals
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: j,l
|
||||
integer,intent(out) :: n_integrals
|
||||
integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num)
|
||||
real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num)
|
||||
|
||||
integer :: i,k
|
||||
double precision :: ao_bielec_integral,cpu_1,cpu_2, wall_1, wall_2
|
||||
double precision :: integral, wall_0
|
||||
double precision :: thresh
|
||||
integer :: kk, m, j1, i1
|
||||
|
||||
thresh = ao_integrals_threshold
|
||||
|
||||
n_integrals = 0
|
||||
|
||||
j1 = j+ishft(l*l-l,-1)
|
||||
do k = 1, ao_num ! r1
|
||||
i1 = ishft(k*k-k,-1)
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
do i = 1, k
|
||||
i1 += 1
|
||||
if (i1 > j1) then
|
||||
exit
|
||||
endif
|
||||
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
|
||||
cycle
|
||||
endif
|
||||
if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then
|
||||
cycle
|
||||
endif
|
||||
!DIR$ FORCEINLINE
|
||||
integral = ao_bielec_integral(i,k,j,l)
|
||||
if (abs(integral) < thresh) then
|
||||
cycle
|
||||
endif
|
||||
n_integrals += 1
|
||||
!DIR$ FORCEINLINE
|
||||
call bielec_integrals_index(i,j,k,l,buffer_i(n_integrals))
|
||||
buffer_value(n_integrals) = integral
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
|
99
src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
Normal file
99
src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
Normal file
@ -0,0 +1,99 @@
|
||||
subroutine ao_bielec_integrals_in_map_slave
|
||||
use map_module
|
||||
use f77_zmq
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes a buffer of integrals
|
||||
END_DOC
|
||||
|
||||
integer :: j,l,n_integrals
|
||||
integer :: rc
|
||||
character*(8), external :: zmq_port
|
||||
integer(ZMQ_PTR) :: zmq_socket_req_inproc, zmq_socket_push_inproc
|
||||
real(integral_kind), allocatable :: buffer_value(:)
|
||||
integer(key_kind), allocatable :: buffer_i(:)
|
||||
|
||||
allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
|
||||
|
||||
! Sockets
|
||||
zmq_socket_req_inproc = f77_zmq_socket(zmq_context, ZMQ_REQ)
|
||||
rc = f77_zmq_connect(zmq_socket_req_inproc, 'inproc://req_rep')
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_socket_req_inproc'
|
||||
endif
|
||||
|
||||
zmq_socket_push_inproc = f77_zmq_socket(zmq_context, ZMQ_PUSH)
|
||||
rc = f77_zmq_connect(zmq_socket_push_inproc, 'inproc://push_pull')
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_socket_push_inproc'
|
||||
endif
|
||||
|
||||
|
||||
|
||||
rc = f77_zmq_send( zmq_socket_req_inproc, 'get_ao_integrals', 16, 0)
|
||||
rc = f77_zmq_recv( zmq_socket_req_inproc, l, 4, 0)
|
||||
|
||||
do while (l > 0)
|
||||
rc = f77_zmq_send( zmq_socket_req_inproc, 'get_ao_integrals', 16, 0)
|
||||
|
||||
do j = 1, l
|
||||
if (ao_overlap_abs(j,l) < ao_integrals_threshold) then
|
||||
cycle
|
||||
endif
|
||||
call compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
|
||||
rc = f77_zmq_send( zmq_socket_push_inproc, n_integrals, 4, ZMQ_SNDMORE)
|
||||
rc = f77_zmq_send( zmq_socket_push_inproc, buffer_i, key_kind*n_integrals, ZMQ_SNDMORE)
|
||||
rc = f77_zmq_send( zmq_socket_push_inproc, buffer_value, integral_kind*n_integrals, 0)
|
||||
enddo
|
||||
rc = f77_zmq_recv( zmq_socket_req_inproc, l, 4, 0)
|
||||
enddo
|
||||
|
||||
deallocate( buffer_i, buffer_value )
|
||||
|
||||
rc = f77_zmq_disconnect(zmq_socket_req_inproc, 'inproc://req_rep')
|
||||
end
|
||||
|
||||
|
||||
subroutine ao_bielec_integrals_in_map_collector
|
||||
use map_module
|
||||
use f77_zmq
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Collects results from the AO integral calculation
|
||||
END_DOC
|
||||
|
||||
integer :: j,l,n_integrals
|
||||
integer :: rc
|
||||
character*(8), external :: zmq_port
|
||||
integer(ZMQ_PTR) :: zmq_socket_pull_inproc
|
||||
real(integral_kind), allocatable :: buffer_value(:)
|
||||
integer(key_kind), allocatable :: buffer_i(:)
|
||||
|
||||
allocate ( buffer_i(ao_num*ao_num), buffer_value(ao_num*ao_num) )
|
||||
|
||||
zmq_socket_pull_inproc = f77_zmq_socket(zmq_context, ZMQ_PULL)
|
||||
rc = f77_zmq_bind(zmq_socket_pull_inproc, 'inproc://push_pull')
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_socket_pull_inproc'
|
||||
endif
|
||||
|
||||
n_integrals = 0
|
||||
do while (n_integrals >= 0)
|
||||
|
||||
rc = f77_zmq_recv( zmq_socket_pull_inproc, n_integrals, 4, 0)
|
||||
if (n_integrals > -1) then
|
||||
rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_i, key_kind*n_integrals, 0)
|
||||
rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_value, integral_kind*n_integrals, 0)
|
||||
call insert_into_ao_integrals_map(n_integrals,buffer_i,buffer_value)
|
||||
else
|
||||
rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_i, key_kind, 0)
|
||||
rc = f77_zmq_recv( zmq_socket_pull_inproc, buffer_value, integral_kind, 0)
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
rc = f77_zmq_unbind(zmq_socket_pull_inproc, 'inproc://push_pull')
|
||||
|
||||
deallocate( buffer_i, buffer_value )
|
||||
end
|
||||
|
@ -247,8 +247,7 @@ BEGIN_PROVIDER [ type(map_type), mo_integrals_map ]
|
||||
print*, 'MO map initialized'
|
||||
END_PROVIDER
|
||||
|
||||
subroutine insert_into_ao_integrals_map(n_integrals, &
|
||||
buffer_i, buffer_values)
|
||||
subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values)
|
||||
use map_module
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -424,15 +423,6 @@ integer*8 function get_mo_map_size()
|
||||
get_mo_map_size = mo_integrals_map % n_elements
|
||||
end
|
||||
|
||||
subroutine clear_mo_map
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Frees the memory of the MO map
|
||||
END_DOC
|
||||
call map_deinit(mo_integrals_map)
|
||||
FREE mo_integrals_map
|
||||
end
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
subroutine dump_$ao_integrals(filename)
|
||||
|
@ -503,3 +503,14 @@ BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_to
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine clear_mo_map
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Frees the memory of the MO map
|
||||
END_DOC
|
||||
call map_deinit(mo_integrals_map)
|
||||
FREE mo_integrals_map mo_bielec_integral_schwartz mo_bielec_integral_jj mo_bielec_integral_jj_anti
|
||||
FREE mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
|
||||
end
|
||||
|
@ -98,7 +98,6 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
|
||||
call write_time(output_mo_basis)
|
||||
|
||||
mo_label = label
|
||||
SOFT_TOUCH mo_coef mo_label
|
||||
end
|
||||
|
||||
subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)
|
||||
|
@ -83,7 +83,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
|
||||
!$OMP DO
|
||||
do i=1,n
|
||||
if ( D(i) < 1.d-12 ) then
|
||||
stop 'Linear dependence in basis set'
|
||||
D(i) = 0.d0
|
||||
else
|
||||
D(i) = 1.d0/dsqrt(D(i))
|
||||
endif
|
||||
|
1
src/ZMQ/NEEDED_CHILDREN_MODULES
Normal file
1
src/ZMQ/NEEDED_CHILDREN_MODULES
Normal file
@ -0,0 +1 @@
|
||||
|
15
src/ZMQ/README.rst
Normal file
15
src/ZMQ/README.rst
Normal file
@ -0,0 +1,15 @@
|
||||
===
|
||||
ZMQ
|
||||
===
|
||||
|
||||
Socket address : defined as an environment variable : QP_RUN_ADDRESS
|
||||
|
||||
|
||||
Needed Modules
|
||||
==============
|
||||
.. Do not edit this section It was auto-generated
|
||||
.. by the `update_README.py` script.
|
||||
Documentation
|
||||
=============
|
||||
.. Do not edit this section It was auto-generated
|
||||
.. by the `update_README.py` script.
|
4
src/ZMQ/f77_zmq_module.f90
Normal file
4
src/ZMQ/f77_zmq_module.f90
Normal file
@ -0,0 +1,4 @@
|
||||
module f77_zmq
|
||||
include 'f77_zmq.h'
|
||||
end module
|
||||
|
105
src/ZMQ/zmq.irp.f
Normal file
105
src/ZMQ/zmq.irp.f
Normal file
@ -0,0 +1,105 @@
|
||||
use f77_zmq
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_context ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Context for the ZeroMQ library
|
||||
END_DOC
|
||||
zmq_context = f77_zmq_ctx_new ()
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ character*(128), qp_run_address ]
|
||||
&BEGIN_PROVIDER [ integer, zmq_port_start ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Address of the qp_run socket
|
||||
! Example : tcp://130.120.229.139:12345
|
||||
END_DOC
|
||||
character*(128) :: buffer
|
||||
call getenv('QP_RUN_ADDRESS',buffer)
|
||||
if (trim(buffer) == '') then
|
||||
stop 'QP_RUN_ADDRESS environment variable not defined'
|
||||
endif
|
||||
|
||||
print *, trim(buffer)
|
||||
integer :: i
|
||||
do i=len(buffer),1,-1
|
||||
if ( buffer(i:i) == ':') then
|
||||
qp_run_address = trim(buffer(1:i-1))
|
||||
read(buffer(i+1:), *) zmq_port_start
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
function zmq_port(ishift)
|
||||
implicit none
|
||||
integer, intent(in) :: ishift
|
||||
character*(8) :: zmq_port
|
||||
write(zmq_port,'(I8)') zmq_port_start+ishift
|
||||
zmq_port = adjustl(trim(zmq_port))
|
||||
end
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_to_qp_run_socket ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Socket on which the qp_run process replies
|
||||
END_DOC
|
||||
integer :: rc
|
||||
zmq_to_qp_run_socket = f77_zmq_socket(zmq_context, ZMQ_REQ)
|
||||
rc = f77_zmq_connect(zmq_to_qp_run_socket, trim(qp_run_address))
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_to_qp_run_socket'
|
||||
endif
|
||||
integer :: i
|
||||
i=4
|
||||
rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_SNDTIMEO, 120000, i)
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to set send timout in zmq_to_qp_run_socket'
|
||||
endif
|
||||
rc = f77_zmq_setsockopt(zmq_to_qp_run_socket, ZMQ_RCVTIMEO, 120000, i)
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to set recv timout in zmq_to_qp_run_socket'
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_push ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Socket on which to push the results (1)
|
||||
END_DOC
|
||||
integer :: rc
|
||||
character*(64) :: address
|
||||
character*(8), external :: zmq_port
|
||||
zmq_socket_push = f77_zmq_socket(zmq_context, ZMQ_PUSH)
|
||||
address = trim(qp_run_address)//':'//zmq_port(1)
|
||||
rc = f77_zmq_connect(zmq_socket_push, trim(address))
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_socket_push'
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(ZMQ_PTR), zmq_socket_pull ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Socket which pulls the results (2)
|
||||
END_DOC
|
||||
integer :: rc
|
||||
character*(64) :: address
|
||||
character*(8), external :: zmq_port
|
||||
zmq_socket_pull = f77_zmq_socket(zmq_context, ZMQ_PULL)
|
||||
address = 'tcp://*:'//zmq_port(2)
|
||||
rc = f77_zmq_bind(zmq_socket_pull, trim(address))
|
||||
if (rc /= 0) then
|
||||
stop 'Unable to connect zmq_socket_pull'
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user