mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-29 20:23:40 +02:00
commit
12f413a389
@ -67,3 +67,15 @@ doc: Use normalized primitive functions
|
||||
interface: ezfio, provider
|
||||
default: true
|
||||
|
||||
[ao_expoim_cosgtos]
|
||||
type: double precision
|
||||
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
|
||||
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
|
||||
interface: ezfio, provider
|
||||
|
||||
[use_cosgtos]
|
||||
type: logical
|
||||
doc: If true, use cosgtos for AO integrals
|
||||
interface: ezfio
|
||||
default: False
|
||||
|
||||
|
33
src/ao_basis/cosgtos.irp.f
Normal file
33
src/ao_basis/cosgtos.irp.f
Normal file
@ -0,0 +1,33 @@
|
||||
BEGIN_PROVIDER [ logical, use_cosgtos ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! If true, use cosgtos for AO integrals
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
if (mpi_master) then
|
||||
call ezfio_has_ao_basis_use_cosgtos(has)
|
||||
if (has) then
|
||||
! write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..'
|
||||
call ezfio_get_ao_basis_use_cosgtos(use_cosgtos)
|
||||
else
|
||||
use_cosgtos = .False.
|
||||
endif
|
||||
endif
|
||||
IRP_IF MPI_DEBUG
|
||||
print *, irp_here, mpi_rank
|
||||
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
|
||||
IRP_ENDIF
|
||||
IRP_IF MPI
|
||||
include 'mpif.h'
|
||||
integer :: ierr
|
||||
call MPI_BCAST( use_cosgtos, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
|
||||
if (ierr /= MPI_SUCCESS) then
|
||||
stop 'Unable to read use_cosgtos with MPI'
|
||||
endif
|
||||
IRP_ENDIF
|
||||
|
||||
! call write_time(6)
|
||||
|
||||
END_PROVIDER
|
@ -1,3 +1,2 @@
|
||||
ao_basis
|
||||
pseudo
|
||||
cosgtos_ao_int
|
||||
|
@ -455,10 +455,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(c)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
ny=0
|
||||
call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in)
|
||||
call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
! call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,R1x,d,nd)
|
||||
else
|
||||
do ix=0,n_pt_in
|
||||
X(ix) = 0.d0
|
||||
@ -469,7 +471,8 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(a-1)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
|
||||
nx = nd
|
||||
do ix=0,n_pt_in
|
||||
@ -479,10 +482,12 @@ recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(c)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
ny=0
|
||||
call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in)
|
||||
call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
! call multiply_poly(Y,ny,R1x,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,R1x,d,nd)
|
||||
endif
|
||||
end
|
||||
|
||||
@ -519,7 +524,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
|
||||
do ix=0,nx
|
||||
X(ix) *= dble(c-1)
|
||||
enddo
|
||||
call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
! call multiply_poly(X,nx,R2x,2,d,nd)
|
||||
call multiply_poly_c2(X,nx,R2x,d,nd)
|
||||
ny = 0
|
||||
do ix=0,dim
|
||||
Y(ix) = 0.d0
|
||||
@ -527,7 +533,8 @@ recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim)
|
||||
|
||||
call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim)
|
||||
if(ny.ge.0)then
|
||||
call multiply_poly(Y,ny,R1xp,2,d,nd)
|
||||
! call multiply_poly(Y,ny,R1xp,2,d,nd)
|
||||
call multiply_poly_c2(Y,ny,R1xp,d,nd)
|
||||
endif
|
||||
endif
|
||||
end
|
||||
|
@ -975,18 +975,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_10,2,d,nd)
|
||||
if (nx >= 0) then
|
||||
integer :: ib
|
||||
do ib=0,nx
|
||||
d(ib ) = d(ib ) + B_10(0) * X(ib)
|
||||
d(ib+1) = d(ib+1) + B_10(1) * X(ib)
|
||||
d(ib+2) = d(ib+2) + B_10(2) * X(ib)
|
||||
enddo
|
||||
|
||||
do nd = nx+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(X,nx,B_10,d,nd)
|
||||
|
||||
nx = nd
|
||||
!DIR$ LOOP COUNT(8)
|
||||
@ -1009,17 +998,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
|
||||
endif
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
if (nx >= 0) then
|
||||
do ib=0,nx
|
||||
d(ib ) = d(ib ) + B_00(0) * X(ib)
|
||||
d(ib+1) = d(ib+1) + B_00(1) * X(ib)
|
||||
d(ib+2) = d(ib+2) + B_00(2) * X(ib)
|
||||
enddo
|
||||
|
||||
do nd = nx+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(X,nx,B_00,d,nd)
|
||||
endif
|
||||
|
||||
ny=0
|
||||
@ -1038,17 +1017,7 @@ recursive subroutine I_x1_pol_mult_recurs(a,c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
if (ny >= 0) then
|
||||
do ib=0,ny
|
||||
d(ib ) = d(ib ) + C_00(0) * Y(ib)
|
||||
d(ib+1) = d(ib+1) + C_00(1) * Y(ib)
|
||||
d(ib+2) = d(ib+2) + C_00(2) * Y(ib)
|
||||
enddo
|
||||
|
||||
do nd = ny+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(Y,ny,C_00,d,nd)
|
||||
end
|
||||
|
||||
recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
@ -1088,18 +1057,7 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
if (nx >= 0) then
|
||||
integer :: ib
|
||||
do ib=0,nx
|
||||
d(ib ) = d(ib ) + B_00(0) * X(ib)
|
||||
d(ib+1) = d(ib+1) + B_00(1) * X(ib)
|
||||
d(ib+2) = d(ib+2) + B_00(2) * X(ib)
|
||||
enddo
|
||||
|
||||
do nd = nx+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(X,nx,B_00,d,nd)
|
||||
|
||||
ny=0
|
||||
|
||||
@ -1111,17 +1069,7 @@ recursive subroutine I_x1_pol_mult_a1(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
if (ny >= 0) then
|
||||
do ib=0,ny
|
||||
d(ib ) = d(ib ) + C_00(0) * Y(ib)
|
||||
d(ib+1) = d(ib+1) + C_00(1) * Y(ib)
|
||||
d(ib+2) = d(ib+2) + C_00(2) * Y(ib)
|
||||
enddo
|
||||
|
||||
do nd = ny+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(Y,ny,C_00,d,nd)
|
||||
|
||||
end
|
||||
|
||||
@ -1150,18 +1098,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_10,2,d,nd)
|
||||
if (nx >= 0) then
|
||||
integer :: ib
|
||||
do ib=0,nx
|
||||
d(ib ) = d(ib ) + B_10(0) * X(ib)
|
||||
d(ib+1) = d(ib+1) + B_10(1) * X(ib)
|
||||
d(ib+2) = d(ib+2) + B_10(2) * X(ib)
|
||||
enddo
|
||||
|
||||
do nd = nx+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(X,nx,B_10,d,nd)
|
||||
|
||||
nx = nd
|
||||
!DIR$ LOOP COUNT(8)
|
||||
@ -1181,17 +1118,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_00,2,d,nd)
|
||||
if (nx >= 0) then
|
||||
do ib=0,nx
|
||||
d(ib ) = d(ib ) + B_00(0) * X(ib)
|
||||
d(ib+1) = d(ib+1) + B_00(1) * X(ib)
|
||||
d(ib+2) = d(ib+2) + B_00(2) * X(ib)
|
||||
enddo
|
||||
|
||||
do nd = nx+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(X,nx,B_00,d,nd)
|
||||
|
||||
ny=0
|
||||
!DIR$ LOOP COUNT(8)
|
||||
@ -1203,17 +1130,7 @@ recursive subroutine I_x1_pol_mult_a2(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_pt_in)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,C_00,2,d,nd)
|
||||
if (ny >= 0) then
|
||||
do ib=0,ny
|
||||
d(ib ) = d(ib ) + C_00(0) * Y(ib)
|
||||
d(ib+1) = d(ib+1) + C_00(1) * Y(ib)
|
||||
d(ib+2) = d(ib+2) + C_00(2) * Y(ib)
|
||||
enddo
|
||||
|
||||
do nd = ny+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(Y,ny,C_00,d,nd)
|
||||
end
|
||||
|
||||
recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
@ -1262,18 +1179,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,D_00,2,d,nd)
|
||||
if (ny >= 0) then
|
||||
integer :: ib
|
||||
do ib=0,ny
|
||||
d(ib ) = d(ib ) + D_00(0) * Y(ib)
|
||||
d(ib+1) = d(ib+1) + D_00(1) * Y(ib)
|
||||
d(ib+2) = d(ib+2) + D_00(2) * Y(ib)
|
||||
enddo
|
||||
|
||||
do nd = ny+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(Y,ny,D_00,d,nd)
|
||||
|
||||
return
|
||||
|
||||
@ -1293,17 +1199,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(X,nx,B_01,2,d,nd)
|
||||
if (nx >= 0) then
|
||||
do ib=0,nx
|
||||
d(ib ) = d(ib ) + B_01(0) * X(ib)
|
||||
d(ib+1) = d(ib+1) + B_01(1) * X(ib)
|
||||
d(ib+2) = d(ib+2) + B_01(2) * X(ib)
|
||||
enddo
|
||||
|
||||
do nd = nx+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(X,nx,B_01,d,nd)
|
||||
|
||||
ny = 0
|
||||
!DIR$ LOOP COUNT(6)
|
||||
@ -1314,17 +1210,7 @@ recursive subroutine I_x2_pol_mult(c,B_10,B_01,B_00,C_00,D_00,d,nd,dim)
|
||||
|
||||
! !DIR$ FORCEINLINE
|
||||
! call multiply_poly(Y,ny,D_00,2,d,nd)
|
||||
if (ny >= 0) then
|
||||
do ib=0,ny
|
||||
d(ib ) = d(ib ) + D_00(0) * Y(ib)
|
||||
d(ib+1) = d(ib+1) + D_00(1) * Y(ib)
|
||||
d(ib+2) = d(ib+2) + D_00(2) * Y(ib)
|
||||
enddo
|
||||
|
||||
do nd = ny+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
endif
|
||||
call multiply_poly_c2(Y,ny,D_00,d,nd)
|
||||
|
||||
end select
|
||||
end
|
||||
|
@ -7,7 +7,13 @@ program bi_ort_ints
|
||||
my_n_pt_r_grid = 10
|
||||
my_n_pt_a_grid = 14
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
call test_3e
|
||||
! call test_3e
|
||||
call test_5idx
|
||||
! call test_5idx2
|
||||
end
|
||||
|
||||
subroutine test_5idx2
|
||||
PROVIDE three_e_5_idx_cycle_2_bi_ort
|
||||
end
|
||||
|
||||
subroutine test_3e
|
||||
@ -16,6 +22,7 @@ subroutine test_3e
|
||||
double precision :: accu, contrib,new,ref
|
||||
i = 1
|
||||
k = 1
|
||||
n = 0
|
||||
accu = 0.d0
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
@ -31,6 +38,7 @@ subroutine test_3e
|
||||
print*,'pb !!'
|
||||
print*,i,k,j,l,m,n
|
||||
print*,ref,new,contrib
|
||||
stop
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
@ -42,3 +50,93 @@ subroutine test_3e
|
||||
|
||||
|
||||
end
|
||||
|
||||
subroutine test_5idx
|
||||
implicit none
|
||||
integer :: i,k,j,l,m,n,ipoint
|
||||
double precision :: accu, contrib,new,ref
|
||||
i = 1
|
||||
k = 1
|
||||
n = 0
|
||||
accu = 0.d0
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
|
||||
new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
|
||||
ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. 1.d-10)then
|
||||
print*,'direct'
|
||||
print*,i,k,j,l,m
|
||||
print*,ref,new,contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
|
||||
ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. 1.d-10)then
|
||||
print*,'exch12'
|
||||
print*,i,k,j,l,m
|
||||
print*,ref,new,contrib
|
||||
stop
|
||||
endif
|
||||
!
|
||||
new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
|
||||
ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. 1.d-10)then
|
||||
print*,'cycle1'
|
||||
print*,i,k,j,l,m
|
||||
print*,ref,new,contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
|
||||
ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. 1.d-10)then
|
||||
print*,'cycle2'
|
||||
print*,i,k,j,l,m
|
||||
print*,ref,new,contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
|
||||
ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. 1.d-10)then
|
||||
print*,'exch23'
|
||||
print*,i,k,j,l,m
|
||||
print*,ref,new,contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
|
||||
ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
|
||||
contrib = dabs(new - ref)
|
||||
accu += contrib
|
||||
if(contrib .gt. 1.d-10)then
|
||||
print*,'exch13'
|
||||
print*,i,k,j,l,m
|
||||
print*,ref,new,contrib
|
||||
stop
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
print*,'accu = ',accu/dble(mo_num)**5
|
||||
|
||||
|
||||
end
|
||||
|
@ -1,7 +1,11 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
@ -14,283 +18,213 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num,
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_direct_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_direct_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
double precision :: wall1, wall0
|
||||
integer :: ipoint
|
||||
double precision, allocatable :: grad_mli(:,:,:), orb_mat(:,:,:)
|
||||
double precision, allocatable :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:,:)
|
||||
double precision, allocatable :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:,:)
|
||||
double precision, allocatable :: tmp_mat(:,:,:,:)
|
||||
allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
|
||||
|
||||
print *, ' Providing the three_e_5_idx_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
do m = 1, mo_num
|
||||
|
||||
allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
|
||||
allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
!$OMP PRIVATE (i,l,ipoint) &
|
||||
!$OMP SHARED (m,mo_num,n_points_final_grid, &
|
||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||
!$OMP grad_mli, orb_mat)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do i=1,mo_num
|
||||
do l=1,mo_num
|
||||
do ipoint=1, n_points_final_grid
|
||||
|
||||
grad_mli(ipoint,l,i) = final_weight_at_r_vector(ipoint) * ( &
|
||||
int2_grad1_u12_bimo_t(ipoint,1,m,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) + &
|
||||
int2_grad1_u12_bimo_t(ipoint,2,m,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) + &
|
||||
int2_grad1_u12_bimo_t(ipoint,3,m,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) )
|
||||
|
||||
orb_mat(ipoint,l,i) = mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, n_points_final_grid, 1.d0, &
|
||||
orb_mat, n_points_final_grid, &
|
||||
grad_mli, n_points_final_grid, 0.d0, &
|
||||
tmp_mat, mo_num*mo_num)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, m, j, i, integral)
|
||||
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = - tmp_mat(l,j,k,i) - tmp_mat(k,i,l,j)
|
||||
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,k,j) - tmp_mat(k,j,l,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_direct_bi_ort', wall1 - wall0
|
||||
deallocate(orb_mat,grad_mli)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = <mlk|-L|jim> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_cycle_1_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_cycle_1_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
|
||||
allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
|
||||
allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num))
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
!$OMP PRIVATE (i,l,ipoint) &
|
||||
!$OMP SHARED (m,mo_num,n_points_final_grid, &
|
||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||
!$OMP rm_grad_ik, lm_grad_ik, rk_grad_im, lk_grad_mi)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do i=1,mo_num
|
||||
do l=1,mo_num
|
||||
do ipoint=1, n_points_final_grid
|
||||
|
||||
lm_grad_ik(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * final_weight_at_r_vector(ipoint)
|
||||
lm_grad_ik(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint)
|
||||
lm_grad_ik(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint)
|
||||
|
||||
rm_grad_ik(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
|
||||
rm_grad_ik(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
|
||||
rm_grad_ik(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
|
||||
|
||||
rk_grad_im(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m)
|
||||
rk_grad_im(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m)
|
||||
rk_grad_im(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
|
||||
lm_grad_ik, 3*n_points_final_grid, &
|
||||
rm_grad_ik, 3*n_points_final_grid, 0.d0, &
|
||||
tmp_mat, mo_num*mo_num)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, j, i, m, integral)
|
||||
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
three_e_5_idx_direct_bi_ort(m,l,j,k,i) = three_e_5_idx_direct_bi_ort(m,l,j,k,i) - tmp_mat(l,j,k,i)
|
||||
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = three_e_5_idx_exch12_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_cycle_1_bi_ort', wall1 - wall0
|
||||
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
|
||||
lm_grad_ik, 3*n_points_final_grid, &
|
||||
rk_grad_im, 3*n_points_final_grid, 0.d0, &
|
||||
tmp_mat, mo_num*mo_num)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = <mlk|-L|imj> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_cycle_2_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_cycle_2_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, i, m, j, integral)
|
||||
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_cycle_2_bi_ort', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_exch23_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_exch23_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, j, m, i, integral)
|
||||
three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = - tmp_mat(l,i,j,k)
|
||||
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = - tmp_mat(k,j,i,l)
|
||||
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = - tmp_mat(k,i,j,l)
|
||||
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = - tmp_mat(l,j,i,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_exch23_bi_ort', wall1 - wall0
|
||||
deallocate(lm_grad_ik)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|ijm> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_exch13_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_exch13_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
!$OMP PRIVATE (i,l,ipoint) &
|
||||
!$OMP SHARED (m,mo_num,n_points_final_grid, &
|
||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||
!$OMP lk_grad_mi)
|
||||
!$OMP DO COLLAPSE(2)
|
||||
do i=1,mo_num
|
||||
do l=1,mo_num
|
||||
do ipoint=1, n_points_final_grid
|
||||
|
||||
lk_grad_mi(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,m,i) * final_weight_at_r_vector(ipoint)
|
||||
lk_grad_mi(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,m,i) * final_weight_at_r_vector(ipoint)
|
||||
lk_grad_mi(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,m,i) * final_weight_at_r_vector(ipoint)
|
||||
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
|
||||
lk_grad_mi, 3*n_points_final_grid, &
|
||||
rm_grad_ik, 3*n_points_final_grid, 0.d0, &
|
||||
tmp_mat, mo_num*mo_num)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, i, j, m, integral)
|
||||
three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(k,j,l,i)
|
||||
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(l,i,k,j)
|
||||
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(l,j,k,i)
|
||||
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(k,i,l,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_exch13_bi_ort', wall1 - wall0
|
||||
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
|
||||
lk_grad_mi, 3*n_points_final_grid, &
|
||||
rk_grad_im, 3*n_points_final_grid, 0.d0, &
|
||||
tmp_mat, mo_num*mo_num)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|mij> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_exch12_bi_ort = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_exch12_bi_ort ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, m, i, j, integral)
|
||||
three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(l,j,i,k)
|
||||
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(k,i,j,l)
|
||||
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(k,j,i,l)
|
||||
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(l,i,j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
deallocate(lk_grad_mi)
|
||||
deallocate(rm_grad_ik)
|
||||
deallocate(rk_grad_im)
|
||||
enddo
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_exch12_bi_ort', wall1 - wall0
|
||||
print *, ' wall time for three_e_5_idx_bi_ort', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
|
295
src/bi_ort_ints/three_body_ijmkl_old.irp.f
Normal file
295
src/bi_ort_ints/three_body_ijmkl_old.irp.f
Normal file
@ -0,0 +1,295 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_direct_bi_ort_old(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_direct_bi_ort_old = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_direct_bi_ort_old ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort_old)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, m, j, i, integral)
|
||||
three_e_5_idx_direct_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_direct_bi_ort_old', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i) = <mlk|-L|jim> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_cycle_1_bi_ort_old = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_cycle_1_bi_ort_old ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort_old)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, j, i, m, integral)
|
||||
three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_cycle_1_bi_ort_old', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i) = <mlk|-L|imj> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_cycle_2_bi_ort_old = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_cycle_2_bi_ort_old ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort_old)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, i, m, j, integral)
|
||||
three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_cycle_2_bi_ort_old', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_exch23_bi_ort_old = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_exch23_bi_ort_old ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort_old)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, j, m, i, integral)
|
||||
three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_exch23_bi_ort_old', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i) = <mlk|-L|ijm> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
three_e_5_idx_exch13_bi_ort_old = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_exch13_bi_ort_old ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort_old)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, i, j, m, integral)
|
||||
three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_exch13_bi_ort_old', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort_old, (mo_num, mo_num, mo_num, mo_num, mo_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
|
||||
!
|
||||
! three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i) = <mlk|-L|mij> ::: notice that i is the RIGHT MO and k is the LEFT MO
|
||||
!
|
||||
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, m, l
|
||||
double precision :: integral, wall1, wall0
|
||||
|
||||
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
|
||||
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
|
||||
|
||||
three_e_5_idx_exch12_bi_ort_old = 0.d0
|
||||
print *, ' Providing the three_e_5_idx_exch12_bi_ort_old ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i,j,k,m,l,integral) &
|
||||
!$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort_old)
|
||||
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
|
||||
do i = 1, mo_num
|
||||
do k = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
do l = 1, mo_num
|
||||
do m = 1, mo_num
|
||||
call give_integrals_3_body_bi_ort(m, l, k, m, i, j, integral)
|
||||
three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i) = -1.d0 * integral
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for three_e_5_idx_exch12_bi_ort_old', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
@ -79,28 +79,31 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
|
||||
integer, intent(in) :: n, l, k, m, j, i
|
||||
double precision, intent(out) :: integral
|
||||
integer :: ipoint
|
||||
double precision :: weight
|
||||
double precision :: weight, tmp
|
||||
|
||||
PROVIDE mo_l_coef mo_r_coef
|
||||
PROVIDE int2_grad1_u12_bimo_t
|
||||
|
||||
integral = 0.d0
|
||||
! (n, l, k, m, j, i)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
weight = final_weight_at_r_vector(ipoint)
|
||||
|
||||
integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
|
||||
tmp = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &
|
||||
* ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j) &
|
||||
+ int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j) &
|
||||
+ int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,l,j) )
|
||||
integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
|
||||
|
||||
tmp = tmp + mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) &
|
||||
* ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i) &
|
||||
+ int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
|
||||
+ int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
|
||||
integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
|
||||
|
||||
tmp = tmp + mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) &
|
||||
* ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i) &
|
||||
+ int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) &
|
||||
+ int2_grad1_u12_bimo_t(ipoint,3,l,j) * int2_grad1_u12_bimo_t(ipoint,3,k,i) )
|
||||
|
||||
integral = integral + tmp * final_weight_at_r_vector(ipoint)
|
||||
enddo
|
||||
|
||||
end subroutine give_integrals_3_body_bi_ort
|
||||
|
@ -198,7 +198,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
|
||||
allocate (bounds(2,nbuckets))
|
||||
do isample=1,nbuckets
|
||||
eta = 1.d0/dble(nbuckets) * dble(isample)
|
||||
ieta = binary_search(waccu,eta,Nabc,ileft,iright)
|
||||
ieta = binary_search(waccu,eta,Nabc)
|
||||
bounds(1,isample) = ileft
|
||||
bounds(2,isample) = ieta
|
||||
ileft = ieta+1
|
||||
|
@ -76,6 +76,8 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset)
|
||||
|
||||
double precision, allocatable :: fock_diag_tmp(:,:)
|
||||
|
||||
if (csubset == 0) return
|
||||
|
||||
allocate(fock_diag_tmp(2,mo_num+1))
|
||||
|
||||
call build_fock_tmp(fock_diag_tmp,psi_det_generators(1,1,i_generator),N_int)
|
||||
@ -177,6 +179,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
|
||||
monoAdo = .true.
|
||||
monoBdo = .true.
|
||||
|
||||
if (csubset == 0) return
|
||||
|
||||
do k=1,N_int
|
||||
hole (k,1) = iand(psi_det_generators(k,1,i_generator), hole_mask(k,1))
|
||||
|
@ -1,19 +0,0 @@
|
||||
[ao_expoim_cosgtos]
|
||||
type: double precision
|
||||
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
|
||||
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
|
||||
interface: ezfio, provider
|
||||
|
||||
[use_cosgtos]
|
||||
type: logical
|
||||
doc: If true, use cosgtos for AO integrals
|
||||
interface: ezfio,provider,ocaml
|
||||
default: False
|
||||
|
||||
[ao_integrals_threshold]
|
||||
type: Threshold
|
||||
doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-15
|
||||
ezfio_name: threshold_ao
|
||||
|
@ -1,2 +0,0 @@
|
||||
ezfio_files
|
||||
ao_basis
|
@ -1,4 +0,0 @@
|
||||
==============
|
||||
cosgtos_ao_int
|
||||
==============
|
||||
|
@ -1,7 +0,0 @@
|
||||
program cosgtos_ao_int
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! TODO : Put the documentation of the program here
|
||||
END_DOC
|
||||
print *, 'Hello world'
|
||||
end
|
@ -6,6 +6,7 @@ BEGIN_PROVIDER [ double precision, cholesky_mo, (mo_num, mo_num, cholesky_ao_num
|
||||
|
||||
integer :: k
|
||||
|
||||
call set_multiple_levels_omp(.False.)
|
||||
print *, 'AO->MO Transformation of Cholesky vectors'
|
||||
!$OMP PARALLEL DO PRIVATE(k)
|
||||
do k=1,cholesky_ao_num
|
||||
|
@ -468,114 +468,6 @@ end subroutine
|
||||
|
||||
|
||||
|
||||
subroutine multiply_poly_0c(b,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:0), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:0+nc)
|
||||
|
||||
integer :: ic
|
||||
|
||||
do ic = 0,nc
|
||||
d(ic) = d(ic) + c(ic) * b(0)
|
||||
enddo
|
||||
|
||||
do nd = nc,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine multiply_poly_1c(b,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:1), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:1+nc)
|
||||
|
||||
integer :: ic, id
|
||||
if(nc < 0) return
|
||||
|
||||
do ic = 0,nc
|
||||
d( ic) = d( ic) + c(ic) * b(0)
|
||||
d(1+ic) = d(1+ic) + c(ic) * b(1)
|
||||
enddo
|
||||
|
||||
do nd = nc+1,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine multiply_poly_2c(b,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:2), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:2+nc)
|
||||
|
||||
integer :: ic, id, k
|
||||
if (nc <0) return
|
||||
|
||||
do ic = 0,nc
|
||||
d( ic) = d( ic) + c(ic) * b(0)
|
||||
d(1+ic) = d(1+ic) + c(ic) * b(1)
|
||||
d(2+ic) = d(2+ic) + c(ic) * b(2)
|
||||
enddo
|
||||
|
||||
do nd = nc+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine multiply_poly_3c(b,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:3), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:3+nc)
|
||||
|
||||
integer :: ic, id
|
||||
if (nc <0) return
|
||||
|
||||
do ic = 1,nc
|
||||
d( ic) = d(1+ic) + c(ic) * b(0)
|
||||
d(1+ic) = d(1+ic) + c(ic) * b(1)
|
||||
d(2+ic) = d(1+ic) + c(ic) * b(2)
|
||||
d(3+ic) = d(1+ic) + c(ic) * b(3)
|
||||
enddo
|
||||
|
||||
do nd = nc+3,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
subroutine multiply_poly(b,nb,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -592,6 +484,30 @@ subroutine multiply_poly(b,nb,c,nc,d,nd)
|
||||
integer :: ib, ic, id, k
|
||||
if(ior(nc,nb) < 0) return !False if nc>=0 and nb>=0
|
||||
|
||||
select case (nb)
|
||||
case (0)
|
||||
call multiply_poly_b0(b,c,nc,d,nd)
|
||||
return
|
||||
case (1)
|
||||
call multiply_poly_b1(b,c,nc,d,nd)
|
||||
return
|
||||
case (2)
|
||||
call multiply_poly_b2(b,c,nc,d,nd)
|
||||
return
|
||||
end select
|
||||
|
||||
select case (nc)
|
||||
case (0)
|
||||
call multiply_poly_c0(b,nb,c,d,nd)
|
||||
return
|
||||
case (1)
|
||||
call multiply_poly_c1(b,nb,c,d,nd)
|
||||
return
|
||||
case (2)
|
||||
call multiply_poly_c2(b,nb,c,d,nd)
|
||||
return
|
||||
end select
|
||||
|
||||
do ib=0,nb
|
||||
do ic = 0,nc
|
||||
d(ib+ic) = d(ib+ic) + c(ic) * b(ib)
|
||||
@ -604,6 +520,254 @@ subroutine multiply_poly(b,nb,c,nc,d,nd)
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine multiply_poly_b0(b,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:0), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:nc)
|
||||
|
||||
integer :: ndtmp
|
||||
integer :: ic, id, k
|
||||
if(nc < 0) return !False if nc>=0
|
||||
|
||||
do ic = 0,nc
|
||||
d(ic) = d(ic) + c(ic) * b(0)
|
||||
enddo
|
||||
|
||||
do nd = nc,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine multiply_poly_b1(b,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:1), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:1+nc)
|
||||
|
||||
integer :: ndtmp
|
||||
integer :: ib, ic, id, k
|
||||
if(nc < 0) return !False if nc>=0
|
||||
|
||||
|
||||
select case (nc)
|
||||
case (0)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1)
|
||||
|
||||
case (1)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
d(2) = d(2) + c(1) * b(1)
|
||||
|
||||
case default
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
do ic = 1,nc
|
||||
d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1)
|
||||
enddo
|
||||
d(nc+1) = d(nc+1) + c(nc) * b(1)
|
||||
|
||||
end select
|
||||
|
||||
do nd = 1+nc,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine multiply_poly_b2(b,c,nc,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nc
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:2), c(0:nc)
|
||||
double precision, intent(inout) :: d(0:2+nc)
|
||||
|
||||
integer :: ndtmp
|
||||
integer :: ib, ic, id, k
|
||||
if(nc < 0) return !False if nc>=0
|
||||
|
||||
select case (nc)
|
||||
case (0)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1)
|
||||
d(2) = d(2) + c(0) * b(2)
|
||||
|
||||
case (1)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
d(2) = d(2) + c(0) * b(2) + c(1) * b(1)
|
||||
d(3) = d(3) + c(1) * b(2)
|
||||
|
||||
case (2)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0)
|
||||
d(3) = d(3) + c(2) * b(1) + c(1) * b(2)
|
||||
d(4) = d(4) + c(2) * b(2)
|
||||
|
||||
case default
|
||||
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
do ic = 2,nc
|
||||
d(ic) = d(ic) + c(ic) * b(0) + c(ic-1) * b(1) + c(ic-2) * b(2)
|
||||
enddo
|
||||
d(nc+1) = d(nc+1) + c(nc) * b(1) + c(nc-1) * b(2)
|
||||
d(nc+2) = d(nc+2) + c(nc) * b(2)
|
||||
|
||||
end select
|
||||
|
||||
do nd = 2+nc,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine multiply_poly_c0(b,nb,c,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nb
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:nb), c(0:0)
|
||||
double precision, intent(inout) :: d(0:nb)
|
||||
|
||||
integer :: ndtmp
|
||||
integer :: ib, ic, id, k
|
||||
if(nb < 0) return !False if nb>=0
|
||||
|
||||
do ib=0,nb
|
||||
d(ib) = d(ib) + c(0) * b(ib)
|
||||
enddo
|
||||
|
||||
do nd = nb,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine multiply_poly_c1(b,nb,c,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nb
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:nb), c(0:1)
|
||||
double precision, intent(inout) :: d(0:nb+1)
|
||||
|
||||
integer :: ndtmp
|
||||
integer :: ib, ic, id, k
|
||||
if(nb < 0) return !False if nb>=0
|
||||
|
||||
select case (nb)
|
||||
case (0)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(1) * b(0)
|
||||
|
||||
case (1)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
d(2) = d(2) + c(1) * b(1)
|
||||
|
||||
case default
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
do ib=1,nb
|
||||
d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1)
|
||||
enddo
|
||||
d(nb+1) = d(nb+1) + c(1) * b(nb)
|
||||
|
||||
end select
|
||||
|
||||
do nd = nb+1,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine multiply_poly_c2(b,nb,c,d,nd)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Multiply two polynomials
|
||||
! D(t) += B(t)*C(t)
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: nb
|
||||
integer, intent(out) :: nd
|
||||
double precision, intent(in) :: b(0:nb), c(0:2)
|
||||
double precision, intent(inout) :: d(0:nb+2)
|
||||
|
||||
integer :: ndtmp
|
||||
integer :: ib, ic, id, k
|
||||
if(nb < 0) return !False if nb>=0
|
||||
|
||||
select case (nb)
|
||||
case (0)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(1) * b(0)
|
||||
d(2) = d(2) + c(2) * b(0)
|
||||
|
||||
case (1)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
d(2) = d(2) + c(1) * b(1) + c(2) * b(0)
|
||||
d(3) = d(3) + c(2) * b(1)
|
||||
|
||||
case (2)
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
d(2) = d(2) + c(0) * b(2) + c(1) * b(1) + c(2) * b(0)
|
||||
d(3) = d(3) + c(1) * b(2) + c(2) * b(1)
|
||||
d(4) = d(4) + c(2) * b(2)
|
||||
|
||||
case default
|
||||
d(0) = d(0) + c(0) * b(0)
|
||||
d(1) = d(1) + c(0) * b(1) + c(1) * b(0)
|
||||
do ib=2,nb
|
||||
d(ib) = d(ib) + c(0) * b(ib) + c(1) * b(ib-1) + c(2) * b(ib-2)
|
||||
enddo
|
||||
d(nb+1) = d(nb+1) + c(1) * b(nb) + c(2) * b(nb-1)
|
||||
d(nb+2) = d(nb+2) + c(2) * b(nb)
|
||||
|
||||
end select
|
||||
|
||||
do nd = nb+2,0,-1
|
||||
if (d(nd) /= 0.d0) exit
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine multiply_poly_v(b,nb,c,nc,d,nd,n_points)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
Loading…
Reference in New Issue
Block a user