mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-09-01 13:43:40 +02:00
beginning to put cholesky in CASSCF
This commit is contained in:
parent
d89682cb7e
commit
d219dc1026
@ -1,7 +1,7 @@
|
|||||||
program spher_harm
|
program spher_harm
|
||||||
implicit none
|
implicit none
|
||||||
! call test_spher_harm
|
call test_spher_harm
|
||||||
! call test_cart
|
! call test_cart
|
||||||
call test_brutal_spheric
|
! call test_brutal_spheric
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -7,6 +7,7 @@ subroutine spher_harm_func_r3(r,l,m,re_ylm, im_ylm)
|
|||||||
double precision :: theta, phi,r_abs
|
double precision :: theta, phi,r_abs
|
||||||
call cartesian_to_spherical(r,theta,phi,r_abs)
|
call cartesian_to_spherical(r,theta,phi,r_abs)
|
||||||
call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
|
call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm)
|
||||||
|
! call spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
@ -131,6 +132,10 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
|
|||||||
tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta)
|
tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta)
|
||||||
re_ylm = tmp * dcos(phi)
|
re_ylm = tmp * dcos(phi)
|
||||||
im_ylm = tmp * dsin(phi)
|
im_ylm = tmp * dsin(phi)
|
||||||
|
else if (l==1.and.m==-1)then
|
||||||
|
tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta)
|
||||||
|
re_ylm = tmp * dcos(phi)
|
||||||
|
im_ylm = -tmp * dsin(phi)
|
||||||
else if(l==1.and.m==0)then
|
else if(l==1.and.m==0)then
|
||||||
tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta)
|
tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta)
|
||||||
re_ylm = tmp
|
re_ylm = tmp
|
||||||
@ -139,10 +144,18 @@ subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm)
|
|||||||
tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
|
tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
|
||||||
re_ylm = tmp * dcos(2.d0*phi)
|
re_ylm = tmp * dcos(2.d0*phi)
|
||||||
im_ylm = tmp * dsin(2.d0*phi)
|
im_ylm = tmp * dsin(2.d0*phi)
|
||||||
|
else if(l==2.and.m==-2)then
|
||||||
|
tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
|
||||||
|
re_ylm = tmp * dcos(2.d0*phi)
|
||||||
|
im_ylm =-tmp * dsin(2.d0*phi)
|
||||||
else if(l==2.and.m==1)then
|
else if(l==2.and.m==1)then
|
||||||
tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
|
tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
|
||||||
re_ylm = tmp * dcos(phi)
|
re_ylm = tmp * dcos(phi)
|
||||||
im_ylm = tmp * dsin(phi)
|
im_ylm = tmp * dsin(phi)
|
||||||
|
else if(l==2.and.m==-1)then
|
||||||
|
tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
|
||||||
|
re_ylm = tmp * dcos(phi)
|
||||||
|
im_ylm =-tmp * dsin(phi)
|
||||||
else if(l==2.and.m==0)then
|
else if(l==2.and.m==0)then
|
||||||
tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0)
|
tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0)
|
||||||
re_ylm = tmp
|
re_ylm = tmp
|
||||||
|
93
src/casscf_cipsi/chol_bielec.irp.f
Normal file
93
src/casscf_cipsi/chol_bielec.irp.f
Normal file
@ -0,0 +1,93 @@
|
|||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, cholesky_no_1_idx_transp, (cholesky_mo_num, n_act_orb, mo_num)]
|
||||||
|
BEGIN_DOC
|
||||||
|
! Cholesky vectors with ONE orbital on the active natural orbital basis
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: i_chol,i_act,i_mo,jj_act
|
||||||
|
double precision, allocatable :: chol_tmp(:,:)
|
||||||
|
allocate(chol_tmp(cholesky_mo_num,n_act_orb))
|
||||||
|
cholesky_no_1_idx_transp = 0.D0
|
||||||
|
do i_mo = 1, mo_num
|
||||||
|
! Get all the integrals corresponding to the "i_mo"
|
||||||
|
do i_act = 1, n_act_orb
|
||||||
|
jj_act = list_act(i_act)
|
||||||
|
do i_chol = 1, cholesky_mo_num
|
||||||
|
chol_tmp(i_chol, i_act) = cholesky_mo_transp(i_chol, jj_act, i_mo)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
! ! Do the matrix product
|
||||||
|
! do i_act = 1, n_act_orb
|
||||||
|
! do jj_act = 1, n_act_orb
|
||||||
|
! do i_chol = 1, cholesky_mo_num
|
||||||
|
! cholesky_no_1_idx_transp(i_chol, i_act, i_mo) += chol_tmp(i_chol, jj_act) * natorbsCI(jj_act,i_act)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0, &
|
||||||
|
chol_tmp, size(chol_tmp,1), &
|
||||||
|
natorbsCI, size(natorbsCI,1), &
|
||||||
|
0.d0, &
|
||||||
|
cholesky_no_1_idx_transp(1,1,i_mo), size(cholesky_no_1_idx_transp,1))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp, (cholesky_mo_num, n_act_orb, n_act_orb)]
|
||||||
|
BEGIN_DOC
|
||||||
|
! Cholesky vectors with TWO orbital on the active natural orbital basis
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: i_chol,i_act,j_act,jj_act
|
||||||
|
double precision, allocatable :: chol_tmp(:,:)
|
||||||
|
allocate(chol_tmp(cholesky_mo_num,n_act_orb))
|
||||||
|
cholesky_no_2_idx_transp = 0.D0
|
||||||
|
do j_act = 1, n_act_orb
|
||||||
|
do i_act = 1, n_act_orb
|
||||||
|
do jj_act = 1, n_act_orb
|
||||||
|
do i_chol = 1, cholesky_mo_num
|
||||||
|
cholesky_no_2_idx_transp(i_chol, i_act, j_act) += cholesky_no_1_idx_transp(i_chol, i_act,jj_act) * natorbsCI(jj_act,i_act)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, cholesky_no_2_idx_transp_dgemm, (cholesky_mo_num, n_act_orb, n_act_orb)]
|
||||||
|
BEGIN_DOC
|
||||||
|
! Cholesky vectors with TWO orbital on the active natural orbital basis
|
||||||
|
END_DOC
|
||||||
|
implicit none
|
||||||
|
integer :: i_chol,i_act,j_act,jj_act
|
||||||
|
double precision, allocatable :: chol_tmp(:,:)
|
||||||
|
allocate(chol_tmp(cholesky_mo_num,n_act_orb))
|
||||||
|
cholesky_no_2_idx_transp_dgemm = 0.D0
|
||||||
|
do j_act = 1, n_act_orb
|
||||||
|
! Get all the integrals corresponding to the "j_act"
|
||||||
|
do i_act = 1, n_act_orb
|
||||||
|
jj_act = list_act(i_act)
|
||||||
|
do i_chol = 1, cholesky_mo_num
|
||||||
|
chol_tmp(i_chol, i_act) = cholesky_no_1_idx_transp(i_chol, j_act, jj_act)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
! ! Do the matrix product
|
||||||
|
! do i_act = 1, n_act_orb
|
||||||
|
! do jj_act = 1, n_act_orb
|
||||||
|
! do i_chol = 1, cholesky_mo_num
|
||||||
|
! cholesky_no_1_idx_transp(i_chol, i_act, j_act) += chol_tmp(i_chol, jj_act) * natorbsCI(jj_act,i_act)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
call dgemm('N','N',cholesky_mo_num,n_act_orb,n_act_orb,1.d0, &
|
||||||
|
chol_tmp, size(chol_tmp,1), &
|
||||||
|
natorbsCI, size(natorbsCI,1), &
|
||||||
|
0.d0, &
|
||||||
|
cholesky_no_2_idx_transp_dgemm(1,1,j_act), size(cholesky_no_2_idx_transp_dgemm,1))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
23
src/casscf_cipsi/test_chol.irp.f
Normal file
23
src/casscf_cipsi/test_chol.irp.f
Normal file
@ -0,0 +1,23 @@
|
|||||||
|
program test_chol
|
||||||
|
implicit none
|
||||||
|
read_wf= .True.
|
||||||
|
touch read_wf
|
||||||
|
call routine
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
subroutine routine
|
||||||
|
implicit none
|
||||||
|
integer :: i_chol, i_act, i_mo
|
||||||
|
double precision :: accu
|
||||||
|
accu = 0.d0
|
||||||
|
do i_mo = 1, n_act_orb
|
||||||
|
do i_act = 1, n_act_orb
|
||||||
|
do i_chol = 1, cholesky_mo_num
|
||||||
|
accu += dabs(cholesky_no_2_idx_transp_dgemm(i_chol,i_act,i_mo) - cholesky_no_2_idx_transp(i_chol,i_act,i_mo))
|
||||||
|
print*,cholesky_no_2_idx_transp_dgemm(i_chol,i_act,i_mo) , cholesky_no_2_idx_transp(i_chol,i_act,i_mo)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
print*,'accu =', accu
|
||||||
|
end
|
@ -101,3 +101,34 @@ BEGIN_PROVIDER [ double precision, cholesky_mo_transp, (cholesky_mo_num, mo_num,
|
|||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, cholesky_semi_mo_transp_simple, (cholesky_mo_num, ao_num, mo_num) ]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Cholesky vectors in MO basis
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
double precision, allocatable :: X(:,:,:)
|
||||||
|
double precision :: wall0, wall1
|
||||||
|
integer :: ierr
|
||||||
|
print *, 'Semi AO->MO Transformation of Cholesky vectors'
|
||||||
|
call wall_time(wall0)
|
||||||
|
|
||||||
|
allocate(X(mo_num,cholesky_mo_num,ao_num), stat=ierr)
|
||||||
|
if (ierr /= 0) then
|
||||||
|
print *, irp_here, ': Allocation failed'
|
||||||
|
endif
|
||||||
|
integer :: i_chol, i_mo, j_mo, i_ao
|
||||||
|
cholesky_semi_mo_transp_simple = 0.d0
|
||||||
|
do i_mo = 1, mo_num
|
||||||
|
do i_ao = 1, ao_num
|
||||||
|
do j_mo = 1, mo_num
|
||||||
|
do i_chol = 1, cholesky_mo_num
|
||||||
|
cholesky_semi_mo_transp_simple(i_chol, i_ao,i_mo) += cholesky_mo_transp(i_chol,j_mo,i_mo) * mo_coef_transp(j_mo,i_ao)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -289,6 +289,106 @@ BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse, (n_points_final_grid)]
|
|||||||
endif
|
endif
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse_bis, (n_points_final_grid)]
|
||||||
|
implicit none
|
||||||
|
integer :: ipoint,m,mm,i,ii,p
|
||||||
|
!!f(R) = \sum_{I} \sum_{J} Phi_I(R) Phi_J(R) V_IJ
|
||||||
|
!! = \sum_{I}\sum_{J}\sum_A Phi_I(R) Phi_J(R) V_AI V_AJ
|
||||||
|
!! = \sum_A \sum_{I}Phi_I(R)V_AI \sum_{J}V_AJ Phi_J(R)
|
||||||
|
!! = \sum_A V_AR G_AR
|
||||||
|
!! V_AR = \sum_{I}Phi_IR V_AI = \sum_{I}Phi^t_RI V_AI
|
||||||
|
double precision :: u_dot_v,wall0,wall1,accu_1, accu_2,mo_i_r1,mo_b_r1
|
||||||
|
double precision :: thresh_1,thresh_2
|
||||||
|
double precision, allocatable :: accu_vec(:),delta_vec(:)
|
||||||
|
thresh_2 = ao_cholesky_threshold * 100.d0
|
||||||
|
thresh_1 = dsqrt(thresh_2)
|
||||||
|
provide cholesky_mo_transp
|
||||||
|
if(elec_alpha_num == elec_beta_num)then
|
||||||
|
call wall_time(wall0)
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE (accu_vec,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) &
|
||||||
|
!$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,mos_in_r_array_omp,aos_in_r_array,thresh_1,thresh_2) &
|
||||||
|
!$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse_bis,n_points_final_grid,cholesky_semi_mo_transp_simple,ao_num)
|
||||||
|
allocate(accu_vec(cholesky_mo_num))
|
||||||
|
!$OMP DO
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
f_hf_cholesky_sparse_bis(ipoint) = 0.d0
|
||||||
|
accu_vec = 0.d0
|
||||||
|
do ii = 1, n_occ_val_orb_for_hf(1)
|
||||||
|
i = list_valence_orb_for_hf(ii,1)
|
||||||
|
mo_i_r1 = mos_in_r_array_omp(i,ipoint)
|
||||||
|
if(dabs(mo_i_r1).lt.thresh_1)cycle
|
||||||
|
do mm = 1, ao_num ! electron 1
|
||||||
|
mo_b_r1 = aos_in_r_array(mm,ipoint)*mo_i_r1
|
||||||
|
if(dabs(mo_b_r1).lt.thresh_2)cycle
|
||||||
|
do p = 1, cholesky_mo_num
|
||||||
|
accu_vec(p) = accu_vec(p) + mo_b_r1 * cholesky_semi_mo_transp_simple(p,mm,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do p = 1, cholesky_mo_num
|
||||||
|
f_hf_cholesky_sparse_bis(ipoint) = f_hf_cholesky_sparse_bis(ipoint) + accu_vec(p) * accu_vec(p)
|
||||||
|
enddo
|
||||||
|
f_hf_cholesky_sparse_bis(ipoint) *= 2.D0
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
deallocate(accu_vec)
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
call wall_time(wall1)
|
||||||
|
print*,'Time to provide f_hf_cholesky_sparse_bis = ',wall1-wall0
|
||||||
|
else
|
||||||
|
call wall_time(wall0)
|
||||||
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
|
!$OMP PRIVATE (accu_vec,delta_vec,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) &
|
||||||
|
!$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,list_basis,mos_in_r_array_omp,thresh_1,thresh_2) &
|
||||||
|
!$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse_bis,n_points_final_grid,cholesky_mo_transp,n_basis_orb)
|
||||||
|
allocate(accu_vec(cholesky_mo_num),delta_vec(cholesky_mo_num))
|
||||||
|
!$OMP DO
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
f_hf_cholesky_sparse_bis(ipoint) = 0.d0
|
||||||
|
accu_vec = 0.d0
|
||||||
|
do ii = 1, n_occ_val_orb_for_hf(2)
|
||||||
|
i = list_valence_orb_for_hf(ii,2)
|
||||||
|
mo_i_r1 = mos_in_r_array_omp(i,ipoint)
|
||||||
|
if(dabs(mo_i_r1).lt.thresh_1)cycle
|
||||||
|
do mm = 1, n_basis_orb ! electron 1
|
||||||
|
m = list_basis(mm)
|
||||||
|
mo_b_r1 = mos_in_r_array_omp(m,ipoint)
|
||||||
|
if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle
|
||||||
|
do p = 1, cholesky_mo_num
|
||||||
|
accu_vec(p) = accu_vec(p) + mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
delta_vec = 0.d0
|
||||||
|
do ii = n_occ_val_orb_for_hf(2)+1,n_occ_val_orb_for_hf(1)
|
||||||
|
i = list_valence_orb_for_hf(ii,1)
|
||||||
|
mo_i_r1 = mos_in_r_array_omp(i,ipoint)
|
||||||
|
if(dabs(mo_i_r1).lt.thresh_1)cycle
|
||||||
|
do mm = 1, n_basis_orb ! electron 1
|
||||||
|
m = list_basis(mm)
|
||||||
|
mo_b_r1 = mos_in_r_array_omp(m,ipoint)
|
||||||
|
if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle
|
||||||
|
do p = 1, cholesky_mo_num
|
||||||
|
delta_vec(p) = delta_vec(p) + mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
do p = 1, cholesky_mo_num
|
||||||
|
f_hf_cholesky_sparse_bis(ipoint) = f_hf_cholesky_sparse_bis(ipoint) + accu_vec(p) * accu_vec(p) + accu_vec(p) * delta_vec(p)
|
||||||
|
enddo
|
||||||
|
f_hf_cholesky_sparse_bis(ipoint) *= 2.D0
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
deallocate(accu_vec)
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
call wall_time(wall1)
|
||||||
|
print*,'Time to provide f_hf_cholesky_sparse_bis = ',wall1-wall0
|
||||||
|
endif
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)]
|
BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: ipoint,i,ii
|
integer :: ipoint,i,ii
|
||||||
|
@ -15,7 +15,23 @@ program projected_operators
|
|||||||
! call test_f_HF_valence_ab
|
! call test_f_HF_valence_ab
|
||||||
! call routine_full_mos
|
! call routine_full_mos
|
||||||
! call test_f_ii_valence_ab
|
! call test_f_ii_valence_ab
|
||||||
call test_f_ia_valence_ab
|
! call test_f_ia_valence_ab
|
||||||
call test_f_ii_ia_aa_valence_ab
|
! call test_f_ii_ia_aa_valence_ab
|
||||||
|
call test
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
subroutine test
|
||||||
|
implicit none
|
||||||
|
integer :: i_point
|
||||||
|
double precision :: ref, new, accu, weight
|
||||||
|
accu = 0.d0
|
||||||
|
do i_point = 1, n_points_final_grid
|
||||||
|
ref = f_hf_cholesky_sparse(i_point)
|
||||||
|
new = f_hf_cholesky_sparse_bis(i_point)
|
||||||
|
weight = final_weight_at_r_vector(i_point)
|
||||||
|
accu += dabs(ref - new) * weight
|
||||||
|
enddo
|
||||||
|
print*,'accu = ',accu
|
||||||
|
|
||||||
|
end
|
||||||
|
Loading…
Reference in New Issue
Block a user