scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-12-22 12:23:37 +01:00
Code Issues Releases Wiki Activity

Curved for Cyrus

Browse Source
This commit is contained in:
Anthony Scemama 2021-06-25 00:03:06 +02:00
parent 08a66a51b5
commit cf3b4fe449
5 changed files with 909 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -6,16 +6,13 @@ program e_curve
integer, allocatable :: iorder(:)
double precision , allocatable :: norm_sort(:)
double precision :: e_0(N_states)
PROVIDE mo_two_e_integrals_in_map
PROVIDE mo_two_e_integrals_in_map mo_one_e_integrals
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
@ -24,19 +21,20 @@ program e_curve
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
PROVIDE psi_bilinear_matrix_values nuclear_repulsion
print *, ''
print *, '=============================='
print *, 'Energies at different cut-offs'
@ -67,27 +65,11 @@ program e_curve
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)
do k=1,N_states
psi_coef(1:N_det,k) = psi_bilinear_matrix_values(1:N_det,k)
call dset_order(psi_coef(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
TOUCH psi_det psi_coef
m = 0
do k=1,n_det
@ -100,10 +82,11 @@ program e_curve
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, &
double precision :: u_dot_u
norm = dsqrt(u_dot_u(psi_coef(1,1),N_det))
print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)', 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
dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, psi_energy(1)
thresh = thresh * dsqrt(10.d0)
enddo
print *, '=========================================================='

63
devel/qmcchem/qmc_e_curve2.irp.f Normal file
View File

@ -0,0 +1,63 @@
program e_curve
use bitmasks
implicit none
integer :: i,j,k, kk, nab, m, l
double precision :: norm, E, hij, num, ci, cj
double precision :: e_0(N_states)
PROVIDE mo_two_e_integrals_in_map mo_one_e_integrals
if (.not.read_wf) then
stop 'Please set read_wf to true'
endif
PROVIDE psi_bilinear_matrix nuclear_repulsion
PROVIDE psi_coef_sorted psi_det psi_coef
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
nab = n_det_alpha_unique+n_det_beta_unique
do while (thresh < 1.d0)
norm = 0.d0
do k=1,n_det
if (dabs(psi_coef(k,1)) < thresh) then
psi_coef(k,1) = 0.d0
endif
norm = norm + psi_coef(k,1)**2
enddo
TOUCH psi_coef
norm = norm/dsqrt(norm)
psi_coef(1:N_det,1) = psi_coef_sorted(1:N_det,1)
psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
do k=1,n_det
if (psi_coef(k,1) == 0.d0) then
exit
endif
enddo
n_det = k-1
TOUCH n_det psi_coef psi_det
j = n_det_alpha_unique+n_det_beta_unique
call u_0_H_u_0(E,psi_coef,n_det,psi_det,N_int,1,size(psi_coef,1))
print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)', thresh, n_det, &
dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-j)), norm, &
psi_energy(1)
thresh = thresh * dsqrt(10.d0)
enddo
print *, '=========================================================='
end

2
devel/trexio/EZFIO.cfg
View File

@ -1,4 +1,4 @@
[trexio_backend]
[backend]
type: integer
doc: Back-end used in TREXIO. 0: HDF5, 1:Text
interface: ezfio, ocaml, provider

23
devel/trexio/export_trexio.irp.f
View File

@ -11,9 +11,9 @@ program export_trexio
print *, 'TREXIO file : '//trim(trexio_filename)
print *, ''
if (trexio_backend == 0) then
if (backend == 0) then
f = trexio_open(trexio_filename, 'w', TREXIO_HDF5)
else if (trexio_backend == 1) then
else if (backend == 1) then
f = trexio_open(trexio_filename, 'w', TREXIO_TEXT)
endif
if (f == 0) then
@ -45,8 +45,8 @@ program export_trexio
rc = trexio_write_nucleus_coord(f, nucl_coord_transp)
call check_success(rc)
! rc = trexio_write_nucleus_label(f, nucl_label)
! call check_success(rc)
rc = trexio_write_nucleus_label(f, nucl_label, 32)
call check_success(rc)
! Pseudo-potentials
@ -90,13 +90,16 @@ program export_trexio
! Basis
! -----
! rc = trexio_write_basis_type(f, 'Gaussian')
! call check_success(rc)
rc = trexio_write_basis_shell_num(f, shell_num)
rc = trexio_write_basis_type(f, 'Gaussian', len('Gaussian'))
call check_success(rc)
rc = trexio_write_basis_shell_center(f, shell_nucl)
rc = trexio_write_basis_num(f, shell_num)
call check_success(rc)
rc = trexio_write_basis_nucleus_shell_num(f, nucleus_shell_num)
call check_success(rc)
rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
call check_success(rc)
rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
@ -119,7 +122,7 @@ program export_trexio
call check_success(rc)
deallocate(factor)
rc = trexio_write_basis_prim_index(f, shell_prim_index)
rc = trexio_write_basis_shell_prim_index(f, shell_prim_index)
call check_success(rc)
rc = trexio_write_basis_exponent(f, prim_expo)

880
devel/trexio/trexio_f.f90
View File

File diff suppressed because it is too large Load Diff