PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-28 11:14:56 +02:00
Code Releases Activity

Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

Browse Source
This commit is contained in:
Anthony Scemama 2025-01-24 11:22:06 +01:00
commit 081df88585
10 changed files with 96 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/becke_numerical_grid/example.irp.f
View File

@ -50,6 +50,7 @@ subroutine example_becke_numerical_grid
print*,'The second example uses the grid points as a collection of spherical grids centered on each atom'
print*,'This is mostly useful if one needs to split contributions between radial/angular/atomic of an integral'
! you browse the nuclei
integral_2 = 0.d0
do i = 1, nucl_num
! you browse the radial points attached to each nucleus
do j = 1, n_points_radial_grid

4
src/becke_numerical_grid/grid_becke_vector.irp.f
View File

@ -71,8 +71,8 @@ END_PROVIDER
enddo
enddo
FREE grid_points_per_atom
FREE final_weight_at_r
! FREE grid_points_per_atom
! FREE final_weight_at_r
call wall_time(wall1)
print *, ' wall time for final_grid_points,', wall1 - wall0

3
src/casscf_cipsi/neworbs.irp.f
View File

@ -116,7 +116,8 @@ END_PROVIDER
end if
end do
if(best_vector_ovrlp_casscf.lt.0)then
best_vector_ovrlp_casscf = minloc(SXeigenval,nMonoEx+1)
! best_vector_ovrlp_casscf = minloc(SXeigenval)
best_vector_ovrlp_casscf = 1
endif
c0=SXeigenvec(1,best_vector_ovrlp_casscf)
if (bavard) then

15
src/determinants/tr_density_matrix.irp.f
View File

@ -1,5 +1,5 @@
BEGIN_PROVIDER [double precision, one_e_tr_dm_mo, (mo_num, mo_num, N_states, N_states)]
use OMP_LIB
implicit none
BEGIN_DOC
@ -16,6 +16,19 @@ BEGIN_PROVIDER [double precision, one_e_tr_dm_mo, (mo_num, mo_num, N_states, N_s
double precision, allocatable :: tmp_a(:,:,:,:), tmp_b(:,:,:,:)
integer :: krow, kcol, lrow, lcol
double precision :: mem_tot_tr_dm
integer :: nthreads
nthreads = OMP_GET_MAX_THREADS()
mem_tot_tr_dm = 8.d0*mo_num*mo_num*n_states*n_states*1d-9*(nthreads+1)
if(mem_tot_tr_dm .gt. 0.9d0 * qp_max_mem) then
print*,'Warning ! you are providing one_e_tr_dm_mo in parallel ! '
print*,'Using that amount of threads ',nthreads
print*,'Each thread needs that amount of Gb ',8.d0*mo_num*mo_num*n_states*n_states*1d-9
print*,'So it takes in total that amount of Gb ',mem_tot_tr_dm
print*,'The max memory for the calculation is ',qp_max_mem
print*,'It may crash because of lack of memory ...'
endif
PROVIDE psi_det_alpha_unique psi_det_beta_unique
one_e_tr_dm_mo = 0d0

23
src/dft_utils_in_r/dipole_becke.irp.f Normal file
View File

@ -0,0 +1,23 @@
BEGIN_PROVIDER [ double precision, mo_dipole_x_becke, (mo_num,mo_num)]
&BEGIN_PROVIDER [ double precision, mo_dipole_y_becke, (mo_num,mo_num)]
&BEGIN_PROVIDER [ double precision, mo_dipole_z_becke, (mo_num,mo_num)]
implicit none
integer :: i,j,ipoint
double precision :: r(3), weight
mo_dipole_x_becke=0.d0
mo_dipole_y_becke=0.d0
mo_dipole_z_becke=0.d0
do i = 1, mo_num
do j = 1, mo_num
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
weight = final_weight_at_r_vector(i)
mo_dipole_x_becke(j,i) += r(1) * weight * mos_in_r_array_transp(ipoint,i) * mos_in_r_array_transp(ipoint,j)
mo_dipole_y_becke(j,i) += r(2) * weight * mos_in_r_array_transp(ipoint,i) * mos_in_r_array_transp(ipoint,j)
mo_dipole_z_becke(j,i) += r(3) * weight * mos_in_r_array_transp(ipoint,i) * mos_in_r_array_transp(ipoint,j)
enddo
enddo
enddo
END_PROVIDER

4
src/iterations/print_summary.irp.f
View File

@ -69,7 +69,9 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
print *, '< S^2 > = ', s2_(k)
print *, 'E = ', e_(k)
print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k)
print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k))
if(dabs(pt2_data % overlap(k,k)).gt.0.d0)then
print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k))
endif
print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
print *, 'rPT2 = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k)
print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)

7
src/mol_properties/EZFIO.cfg
View File

@ -1,3 +1,10 @@
[dft_grid_mo_dipole]
type: logical
doc: If true, you compute the MO x/y/z dipole with DFT's grid
interface: ezfio,provider,ocaml
default: false
[print_all_transitions]
type: logical
doc: If true, print the transition between all the states

1
src/mol_properties/NEED
View File

@ -1,2 +1,3 @@
determinants
davidson_undressed
dft_utils_in_r

14
src/mol_properties/mo_dipole_properties.irp.f Normal file
View File

@ -0,0 +1,14 @@
BEGIN_PROVIDER [double precision, mo_prop_dipole_x , (mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, mo_prop_dipole_y , (mo_num,mo_num)]
&BEGIN_PROVIDER [double precision, mo_prop_dipole_z , (mo_num,mo_num)]
implicit none
if(dft_grid_mo_dipole)then
mo_prop_dipole_x=mo_dipole_x_becke
mo_prop_dipole_y=mo_dipole_y_becke
mo_prop_dipole_z=mo_dipole_z_becke
else
mo_prop_dipole_x=mo_dipole_x
mo_prop_dipole_y=mo_dipole_y
mo_prop_dipole_z=mo_dipole_z
endif
END_PROVIDER

53
src/mol_properties/multi_s_dipole_moment.irp.f
View File

@ -39,24 +39,28 @@
! p,q: general spatial MOs
! gamma^{nm}: density matrix \bra{\Psi^n} a^{\dagger}_a a_i \ket{\Psi^m}
END_DOC
USE OMP_LIB
integer :: istate, jstate ! States
integer :: i, j ! general spatial MOs
double precision :: nuclei_part_x, nuclei_part_y, nuclei_part_z
double precision :: mem_tot_tr_dm
integer :: nthreads
nthreads = OMP_GET_MAX_THREADS()
mem_tot_tr_dm = 8.d0*mo_num*mo_num*n_states*n_states*1d-9*(nthreads+1)
multi_s_x_dipole_moment = 0.d0
multi_s_y_dipole_moment = 0.d0
multi_s_z_dipole_moment = 0.d0
if(8.d0*mo_num*mo_num*n_states*n_states*1d-9 .lt. 200.d0) then
if(mem_tot_tr_dm .lt. 0.9d0 * qp_max_mem) then
do jstate = 1, N_states
do istate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
multi_s_x_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_x(j,i)
multi_s_y_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_y(j,i)
multi_s_z_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_dipole_z(j,i)
multi_s_x_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_prop_dipole_x(j,i)
multi_s_y_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_prop_dipole_y(j,i)
multi_s_z_dipole_moment(istate,jstate) -= one_e_tr_dm_mo(j,i,istate,jstate) * mo_prop_dipole_z(j,i)
enddo
enddo
enddo
@ -66,6 +70,7 @@
! no enouph memory
! on the fly scheme
print*,'Computing on the fly the various dipole matrix elements '
PROVIDE psi_det_alpha_unique psi_det_beta_unique
@ -85,7 +90,7 @@
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique, &
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values, &
!$OMP mo_dipole_x, mo_dipole_y, mo_dipole_z, &
!$OMP mo_prop_dipole_x, mo_prop_dipole_y, mo_prop_dipole_z, &
!$OMP multi_s_x_dipole_moment, multi_s_y_dipole_moment, multi_s_z_dipole_moment)
!$OMP DO COLLAPSE(2)
do istate = 1, N_states
@ -103,9 +108,9 @@
ck = psi_bilinear_matrix_values(k_a,istate)*psi_bilinear_matrix_values(k_a,jstate)
do l = 1, elec_alpha_num
j = occ(l,1)
multi_s_x_dipole_moment(istate,jstate) -= ck * mo_dipole_x(j,j)
multi_s_y_dipole_moment(istate,jstate) -= ck * mo_dipole_y(j,j)
multi_s_z_dipole_moment(istate,jstate) -= ck * mo_dipole_z(j,j)
multi_s_x_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_x(j,j)
multi_s_y_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_y(j,j)
multi_s_z_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_z(j,j)
enddo
if (k_a == N_det) cycle
@ -121,13 +126,13 @@
call get_single_excitation_spin(tmp_det(1,1), tmp_det2, exc, phase, N_int)
call decode_exc_spin(exc, h1, p1, h2, p2)
ckl = psi_bilinear_matrix_values(k_a,istate)*psi_bilinear_matrix_values(l,jstate) * phase
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(h1,p1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(h1,p1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(h1,p1)
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(h1,p1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(h1,p1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(h1,p1)
ckl = psi_bilinear_matrix_values(k_a,jstate)*psi_bilinear_matrix_values(l,istate) * phase
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(p1,h1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(p1,h1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(p1,h1)
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(p1,h1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(p1,h1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(p1,h1)
endif
l = l+1
if (l > N_det) exit
@ -148,9 +153,9 @@
ck = psi_bilinear_matrix_transp_values(k_b,istate)*psi_bilinear_matrix_transp_values(k_b,jstate)
do l = 1, elec_beta_num
j = occ(l,2)
multi_s_x_dipole_moment(istate,jstate) -= ck * mo_dipole_x(j,j)
multi_s_y_dipole_moment(istate,jstate) -= ck * mo_dipole_y(j,j)
multi_s_z_dipole_moment(istate,jstate) -= ck * mo_dipole_z(j,j)
multi_s_x_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_x(j,j)
multi_s_y_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_y(j,j)
multi_s_z_dipole_moment(istate,jstate) -= ck * mo_prop_dipole_z(j,j)
enddo
if (k_b == N_det) cycle
@ -166,13 +171,13 @@
call get_single_excitation_spin(tmp_det(1,2), tmp_det2, exc, phase, N_int)
call decode_exc_spin(exc, h1, p1, h2, p2)
ckl = psi_bilinear_matrix_transp_values(k_b,istate)*psi_bilinear_matrix_transp_values(l,jstate) * phase
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(h1,p1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(h1,p1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(h1,p1)
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(h1,p1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(h1,p1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(h1,p1)
ckl = psi_bilinear_matrix_transp_values(k_b,jstate)*psi_bilinear_matrix_transp_values(l,istate) * phase
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_dipole_x(p1,h1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_dipole_y(p1,h1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_dipole_z(p1,h1)
multi_s_x_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_x(p1,h1)
multi_s_y_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_y(p1,h1)
multi_s_z_dipole_moment(istate,jstate) -= ckl * mo_prop_dipole_z(p1,h1)
endif
l = l+1
if (l > N_det) exit