scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2025-01-09 04:43:07 +01:00
Code Issues Releases Wiki Activity

Compare commits

base: scemama:dd65b2677b1a88d3ad2cfa34bf3362237b879919
Branches Tags
scemama:master
..
compare: scemama:5e60f68b89954071e4885d76a681de9cff2328e3
Branches Tags
scemama:master

No commits in common. "dd65b2677b1a88d3ad2cfa34bf3362237b879919" and "5e60f68b89954071e4885d76a681de9cff2328e3" have entirely different histories.
dd65b2677b ... 5e60f68b89

16 changed files with 0 additions and 593 deletions
Show Stats Expand all files Collapse all files

25
devel/cassdtq/EZFIO.cfg
View File

@ -1,25 +0,0 @@
[energy]
type: double precision
doc: Calculated Selected CASSD energy
interface: ezfio
size: (determinants.n_states)
[energy_pt2]
type: double precision
doc: Calculated CASSD energy + PT2
interface: ezfio
size: (determinants.n_states)
[do_ddci]
type: logical
doc: If true, remove purely inactive double excitations
interface: ezfio,provider,ocaml
default: False
[do_only_1h1p]
type: logical
doc: If true, do only one hole/one particle excitations
interface: ezfio,provider,ocaml
default: False

3
devel/cassdtq/NEED
View File

@ -1,3 +0,0 @@
cipsi
selectors_full

5
devel/cassdtq/README.rst
View File

@ -1,5 +0,0 @@
======
cisdtq
======
Selected CI in the CISDTQ space.

22
devel/cassdtq/cassdtq.irp.f
View File

@ -1,22 +0,0 @@
program cassdtq
implicit none
BEGIN_DOC
! Selected CISDTQ with stochastic selection and PT2.
END_DOC
if (.not.is_zmq_slave) then
PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
if (do_pt2) then
call run_stochastic_cipsi
else
call run_cipsi
endif
else
PROVIDE mo_two_e_integrals_in_map
call run_slave_cipsi
endif
end

8
devel/cassdtq/class.irp.f
View File

@ -1,8 +0,0 @@
BEGIN_PROVIDER [ logical, do_only_cas ]
implicit none
BEGIN_DOC
! In the CASSDTQ case, all those are always false
END_DOC
do_only_cas = .False.
END_PROVIDER

78
devel/cassdtq/generators.irp.f
View File

@ -1,78 +0,0 @@
use bitmasks
logical function is_a_generator(det)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: det(N_int,2)
integer :: na, nb
integer, external :: number_of_holes, number_of_particles
provide N_int
na = number_of_holes(det)
nb = number_of_particles(det)
is_a_generator = (na <= 2) .and. (nb <= 2)
end
BEGIN_PROVIDER [ integer, N_det_generators ]
implicit none
BEGIN_DOC
! Number of generator detetrminants
END_DOC
integer :: i,k,l
logical, external :: is_a_generator
provide N_int
call write_time(6)
N_det_generators = 0
do i=1,N_det
if (is_a_generator(psi_det_sorted(1,1,i))) then
N_det_generators += 1
endif
enddo
N_det_generators = max(N_det_generators,1)
call write_int(6,N_det_generators,'Number of generators')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
&BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_gen, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_sorted_gen, (psi_det_size,N_states) ]
&BEGIN_PROVIDER [ integer, psi_det_sorted_gen_order, (psi_det_size) ]
implicit none
BEGIN_DOC
! For Single reference wave functions, the generator is the
! Hartree-Fock determinant
END_DOC
integer :: i, k, l, m
logical, external :: is_a_generator
integer, allocatable :: nongen(:)
integer :: inongen
allocate(nongen(N_det))
inongen = 0
m=0
do i=1,N_det
if (is_a_generator(psi_det_sorted(1,1,i))) then
m = m+1
psi_det_sorted_gen_order(i) = m
do k=1,N_int
psi_det_generators(k,1,m) = psi_det_sorted(k,1,i)
psi_det_generators(k,2,m) = psi_det_sorted(k,2,i)
enddo
psi_coef_generators(m,:) = psi_coef_sorted(i,:)
else
inongen += 1
nongen(inongen) = i
endif
enddo
psi_det_sorted_gen(:,:,:N_det_generators) = psi_det_generators(:,:,:N_det_generators)
psi_coef_sorted_gen(:N_det_generators, :) = psi_coef_generators(:N_det_generators, :)
do i=1,inongen
psi_det_sorted_gen_order(nongen(i)) = N_det_generators+i
psi_det_sorted_gen(:,:,N_det_generators+i) = psi_det_sorted(:,:,nongen(i))
psi_coef_sorted_gen(N_det_generators+i, :) = psi_coef_sorted(nongen(i),:)
end do
END_PROVIDER

9
devel/cassdtq/save_energy.irp.f
View File

@ -1,9 +0,0 @@
subroutine save_energy(E,pt2)
implicit none
BEGIN_DOC
! Saves the energy in |EZFIO|.
END_DOC
double precision, intent(in) :: E(N_states), pt2(N_states)
call ezfio_set_cassdtq_energy(E(1:N_states))
call ezfio_set_cassdtq_energy_pt2(E(1:N_states)+pt2(1:N_states))
end

6
devel/qmcchem/EZFIO.cfg
View File

@ -1,6 +0,0 @@
[ci_threshold]
type: Threshold
doc: Threshold on the CI coefficients as computed in QMCChem
interface: ezfio,provider,ocaml
default: 1.e-8

4
devel/qmcchem/NEED
View File

@ -1,4 +0,0 @@
determinants
davidson_undressed
fci
zmq

12
devel/qmcchem/README.rst
View File

@ -1,12 +0,0 @@
==============
QmcChem Module
==============
For multi-state calculations, to extract state 2 use:
``
QP_STATE=2 qp_run save_for_qmcchem x.ezfio
``
(state 1 is the ground state).

8
devel/qmcchem/expand_spindets_qmcchem.irp.f
View File

@ -1,8 +0,0 @@
program densify
implicit none
read_wf = .True.
touch read_wf
call generate_all_alpha_beta_det_products()
call diagonalize_ci
call save_wavefunction
end

131
devel/qmcchem/pot_ao_pseudo_ints.irp.f
View File

@ -1,131 +0,0 @@
BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
!
! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
END_DOC
! l,m : Y(l,m) parameters
! c(3) : pseudopotential center
! a(3) : Atomic Orbital center
! n_a(3) : Powers of x,y,z in the Atomic Orbital
! g_a : Atomic Orbital exponent
! r : Distance between the Atomic Orbital center and the considered point
double precision, external :: ylm_orb
integer :: n_a(3)
double precision :: a(3), c(3), g_a
integer :: i,j,k,l,m,n,p
double precision :: dr, Ulc
double precision :: y
double precision, allocatable :: r(:)
allocate (r(pseudo_grid_size))
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
do j=2,pseudo_grid_size
r(j) = r(j-1) + dr
enddo
ao_pseudo_grid = 0.d0
do j=1,pseudo_grid_size
do k=1,nucl_num
c(1:3) = nucl_coord(k,1:3)
do l=0,pseudo_lmax
do i=1,ao_num
a(1:3) = nucl_coord(ao_nucl(i),1:3)
n_a(1:3) = ao_power(i,1:3)
do p=1,ao_prim_num(i)
g_a = ao_expo_ordered_transp(p,i)
do m=-l,l
y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
ao_pseudo_grid(i,m,l,k,j) = ao_pseudo_grid(i,m,l,k,j) + &
ao_coef_normalized_ordered_transp(p,i)*y
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(r)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
!
! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
END_DOC
! l,m : Y(l,m) parameters
! c(3) : pseudopotential center
! a(3) : Atomic Orbital center
! n_a(3) : Powers of x,y,z in the Atomic Orbital
! g_a : Atomic Orbital exponent
! r : Distance between the Atomic Orbital center and the considered point
double precision, external :: ylm_orb
integer :: n_a(3)
double precision :: a(3), c(3), g_a
integer :: i,j,k,l,m,n,p
double precision :: dr, Ulc
double precision :: y
double precision, allocatable :: r(:)
allocate (r(pseudo_grid_size))
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
do j=2,pseudo_grid_size
r(j) = r(j-1) + dr
enddo
mo_pseudo_grid = 0.d0
do n=1,pseudo_grid_size
do k=1,nucl_num
do l=0,pseudo_lmax
do m=-l,l
do i=1,ao_num
do j=1,mo_num
if (dabs(ao_pseudo_grid(i,m,l,k,n)) < 1.e-12) then
cycle
endif
if (dabs(mo_coef(i,j)) < 1.e-8) then
cycle
endif
mo_pseudo_grid(j,m,l,k,n) = mo_pseudo_grid(j,m,l,k,n) + &
ao_pseudo_grid(i,m,l,k,n) * mo_coef(i,j)
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(r)
END_PROVIDER
double precision function test_pseudo_grid_ao(i,j)
implicit none
integer, intent(in) :: i,j
integer :: k,l,m,n
double precision :: r, dr,u
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
test_pseudo_grid_ao = 0.d0
r = 0.d0
do k=1,pseudo_grid_size
do n=1,nucl_num
do l = 0,pseudo_lmax
u = pseudo_v_kl(n,l,1) * exp(-pseudo_dz_kl(n,l,1)*r*r)* r*r*dr
do m=-l,l
test_pseudo_grid_ao += ao_pseudo_grid(i,m,l,n,k) * ao_pseudo_grid(j,m,l,n,k) * u
enddo
enddo
enddo
r = r+dr
enddo
end

8
devel/qmcchem/pseudo.irp.f
View File

@ -1,8 +0,0 @@
subroutine write_pseudopotential
implicit none
BEGIN_DOC
! Write the pseudo_potential into the EZFIO file
END_DOC
call ezfio_set_pseudo_mo_pseudo_grid(mo_pseudo_grid)
end

113
devel/qmcchem/qmc_e_curve.irp.f
View File

@ -1,113 +0,0 @@
program e_curve
use bitmasks
implicit none
integer :: i,j,k, kk, nab, m, l
double precision :: norm, E, hij, num, ci, cj
integer, allocatable :: iorder(:)
double precision , allocatable :: norm_sort(:)
double precision :: e_0(N_states)
PROVIDE mo_two_e_integrals_in_map
nab = n_det_alpha_unique+n_det_beta_unique
allocate ( norm_sort(0:nab), iorder(0:nab) )
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
double precision, allocatable :: u_0(:,:), v_0(:,:)
allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
allocate(u_0(N_states,N_det),v_0(N_states,N_det))
norm_sort(0) = 0.d0
iorder(0) = 0
do i=1,n_det_alpha_unique
norm_sort(i) = det_alpha_norm(i)
iorder(i) = i
enddo
do i=1,n_det_beta_unique
norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
iorder(i+n_det_alpha_unique) = -i
enddo
call dsort(norm_sort(1),iorder(1),nab)
if (.not.read_wf) then
stop 'Please set read_wf to true'
endif
PROVIDE psi_bilinear_matrix_values nuclear_repulsion
print *, ''
print *, '=============================='
print *, 'Energies at different cut-offs'
print *, '=============================='
print *, ''
print *, '=========================================================='
print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)', 'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
print *, '=========================================================='
double precision :: thresh
integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
thresh = 1.d-10
do j=0,nab
i = iorder(j)
if (i<0) then
do k=1,n_det
if (psi_bilinear_matrix_columns(k) == -i) then
psi_bilinear_matrix_values(k,1) = 0.d0
endif
enddo
else
do k=1,n_det
if (psi_bilinear_matrix_rows(k) == i) then
psi_bilinear_matrix_values(k,1) = 0.d0
endif
enddo
endif
if (thresh > norm_sort(j)) then
cycle
endif
u_0 = psi_bilinear_matrix_values(1:N_det,1:N_states)
v_t = 0.d0
s_t = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_states)
call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_states, N_det)
double precision, external :: u_dot_u, u_dot_v
do i=1,N_states
e_0(i) = u_dot_v(v_t(1,i),u_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
enddo
m = 0
do k=1,n_det
if (psi_bilinear_matrix_values(k,1) /= 0.d0) then
m = m+1
endif
enddo
if (m == 0) then
exit
endif
E = E_0(1) + nuclear_repulsion
norm = u_dot_u(u_0(1,1),N_det)
print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)', thresh, m, &
dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, E
thresh = thresh * dsqrt(10.d0)
enddo
print *, '=========================================================='
deallocate (iorder, norm_sort)
end

45
devel/qmcchem/save_for_qmcchem.irp.f
View File

@ -1,45 +0,0 @@
program save_for_qmc
integer :: iunit
logical :: exists
double precision :: e_ref
! Determinants
read_wf = .True.
TOUCH read_wf
print *, "N_det = ", N_det
call write_spindeterminants
! Reference Energy
if (do_pseudo) then
call write_pseudopotential
endif
call system( &
'mkdir -p '//trim(ezfio_filename)//'/simulation ;' // &
'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/simulation/.version ; ' // &
'mkdir -p '//trim(ezfio_filename)//'/properties ;' // &
'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/properties/.version ; ' // &
'echo T > '//trim(ezfio_filename)//'/properties/e_loc' &
)
iunit = 13
open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write')
call ezfio_has_fci_energy_pt2(exists)
if (exists) then
call ezfio_get_fci_energy_pt2(e_ref)
else
call ezfio_has_fci_energy(exists)
if (exists) then
call ezfio_get_fci_energy(e_ref)
else
call ezfio_has_hartree_fock_energy(exists)
if (exists) then
call ezfio_get_hartree_fock_energy(e_ref)
else
e_ref = 0.d0
endif
endif
endif
write(iunit,*) e_ref
close(iunit)
end

116
devel/qmcchem/truncate_wf_qmcchem.irp.f
View File

@ -1,116 +0,0 @@
program truncate
read_wf = .True.
SOFT_TOUCH read_wf
call run
end
subroutine run
use bitmasks
implicit none
integer :: i,j,k, kk, nab, m, l
double precision :: norm, E, hij, num, ci, cj
integer, allocatable :: iorder(:)
double precision , allocatable :: norm_sort(:)
double precision :: e_0(N_states)
PROVIDE mo_two_e_integrals_in_map H_apply_buffer_allocated
nab = n_det_alpha_unique+n_det_beta_unique
allocate ( norm_sort(0:nab), iorder(0:nab) )
integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
double precision, allocatable :: u_t(:,:), v_t(:,:), s_t(:,:)
double precision, allocatable :: u_0(:,:), v_0(:,:)
allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
allocate(u_0(N_det,N_states),v_0(N_det,N_states))
norm_sort(0) = 0.d0
iorder(0) = 0
do i=1,n_det_alpha_unique
norm_sort(i) = det_alpha_norm(i)
iorder(i) = i
enddo
do i=1,n_det_beta_unique
norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
iorder(i+n_det_alpha_unique) = -i
enddo
call dsort(norm_sort(1),iorder(1),nab)
PROVIDE psi_bilinear_matrix_values psi_bilinear_matrix_rows psi_bilinear_matrix_columns
PROVIDE nuclear_repulsion
print *, ''
do j=0,nab
i = iorder(j)
if (i<0) then
!$OMP PARALLEL DO PRIVATE(k)
do k=1,n_det
if (psi_bilinear_matrix_columns(k) == -i) then
do l=1,N_states
psi_bilinear_matrix_values(k,l) = 0.d0
enddo
endif
enddo
!$OMP END PARALLEL DO
else
!$OMP PARALLEL DO PRIVATE(k)
do k=1,n_det
if (psi_bilinear_matrix_rows(k) == i) then
do l=1,N_states
psi_bilinear_matrix_values(k,l) = 0.d0
enddo
endif
enddo
!$OMP END PARALLEL DO
endif
if (ci_threshold <= norm_sort(j)) then
exit
endif
enddo
m = 0
do k=1,n_det
if (sum(psi_bilinear_matrix_values(k,1:N_states)) /= 0.d0) then
m = m+1
endif
enddo
do k=1,N_states
E = E_0(k) + nuclear_repulsion
enddo
print *, 'Number of determinants:', m
call wf_of_psi_bilinear_matrix(.True.)
call save_wavefunction
u_0(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
v_0(1:N_det,1:N_states) = 0.d0
u_t(1:N_states,1:N_det) = 0.d0
v_t(1:N_states,1:N_det) = 0.d0
s_t(1:N_states,1:N_det) = 0.d0
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_states)
print *, 'Computing H|Psi> ...'
call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
print *, 'Done'
call dtranspose( &
v_t, &
size(v_t, 1), &
v_0, &
size(v_0, 1), &
N_states, N_det)
double precision, external :: u_dot_u, u_dot_v
do i=1,N_states
e_0(i) = u_dot_v(u_0(1,i),v_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
print *, 'E(',i,') = ', e_0(i) + nuclear_repulsion
enddo
deallocate (iorder, norm_sort)
end