LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-03 20:54:00 +01:00
Code Issues Releases Wiki Activity

Work on QMC=Chem interface

Browse Source
This commit is contained in:
Anthony Scemama 2015-06-21 23:08:18 +02:00
parent f15c6583ec
commit 9548cf14d4
4 changed files with 6 additions and 389 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/gfortran.cfg
View File

@ -11,7 +11,7 @@
#
[COMMON]
FC : gfortran -ffree-line-length-none -I .
LAPACK_LIB : -lblas -llapack
LAPACK_LIB : -llapack -lblas
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32

2
scripts/compilation/qp_create_ninja.py
View File

@ -249,7 +249,7 @@ def ninja_ezfio_rule():
install_lib_ezfio = join(QP_ROOT, 'install', 'EZFIO', "lib", "libezfio.a")
l_cmd = ["cd {0}".format(QP_ROOT_EZFIO)] + l_flag
l_cmd += ["ninja && ln -f {0} {1}".format(install_lib_ezfio, EZFIO_LIB)]
l_cmd += ["ninja && ln -sf {0} {1}".format(install_lib_ezfio, EZFIO_LIB)]
l_string = ["rule build_ezfio",
" command = {0}".format(" ; ".join(l_cmd)),

228
src/QmcChem/pot_ao_pseudo_ints.irp.f
View File

@ -1,226 +1,3 @@
BEGIN_PROVIDER [ double precision, aux_pseudo_integral, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
implicit none
BEGIN_DOC
! Pseudo-potential
END_DOC
if (do_pseudo) then
! aux_pseudo_integral = aux_pseudo_integral_local + aux_pseudo_integral_non_local
! aux_pseudo_integral = aux_pseudo_integral_local
aux_pseudo_integral = aux_pseudo_integral_non_local
else
aux_pseudo_integral = 0.d0
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, aux_pseudo_integral_local, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
implicit none
BEGIN_DOC
! Local pseudo-potential
END_DOC
double precision :: alpha, beta, gama, delta
integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt_in,m
double precision :: Vloc, Vpseudo
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
integer :: thread_num
aux_pseudo_integral_local = 0.d0
!! Dump array
integer, allocatable :: n_k_dump(:)
double precision, allocatable :: v_k_dump(:), dz_k_dump(:)
allocate(n_k_dump(1:pseudo_klocmax), v_k_dump(1:pseudo_klocmax), dz_k_dump(1:pseudo_klocmax))
! _
! / _. | _ |
! \_ (_| | (_ |_| |
!
print*, 'Providing the nuclear electron pseudo integrals '
call wall_time(wall_1)
call cpu_time(cpu_1)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
!$OMP num_A,num_B,Z,c,n_pt_in, &
!$OMP v_k_dump,n_k_dump, dz_k_dump, &
!$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (aux_basis_num_sqrt,aux_basis_prim_num,aux_basis_expo_transp,aux_basis_power,aux_basis_nucl,nucl_coord,aux_basis_coef_transp, &
!$OMP aux_pseudo_integral_local,nucl_num,nucl_charge, &
!$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k,pseudo_n_k, pseudo_dz_k, &
!$OMP wall_1)
!$OMP DO SCHEDULE (guided)
do j = 1, aux_basis_num_sqrt
num_A = aux_basis_nucl(j)
power_A(1:3)= aux_basis_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, aux_basis_num_sqrt
num_B = aux_basis_nucl(i)
power_B(1:3)= aux_basis_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,aux_basis_prim_num(j)
alpha = aux_basis_expo_transp(l,j)
do m=1,aux_basis_prim_num(i)
beta = aux_basis_expo_transp(m,i)
double precision :: c
c = 0.d0
do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
v_k_dump = pseudo_v_k(k,1:pseudo_klocmax)
n_k_dump = pseudo_n_k(k,1:pseudo_klocmax)
dz_k_dump = pseudo_dz_k(k,1:pseudo_klocmax)
c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump, &
A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
aux_pseudo_integral_local(i,j) = aux_pseudo_integral_local(i,j) + &
aux_basis_coef_transp(l,j)*aux_basis_coef_transp(m,i)*c
enddo
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(j)/float(aux_basis_num_sqrt), '% in ', &
wall_2-wall_1, 's'
endif
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
deallocate(n_k_dump,v_k_dump, dz_k_dump)
END_PROVIDER
BEGIN_PROVIDER [ double precision, aux_pseudo_integral_non_local, (aux_basis_num_sqrt,aux_basis_num_sqrt)]
implicit none
BEGIN_DOC
! Local pseudo-potential
END_DOC
double precision :: alpha, beta, gama, delta
integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(3)
integer :: i,j,k,l,n_pt_in,m
double precision :: Vloc, Vpseudo
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
integer :: thread_num
aux_pseudo_integral_non_local = 0.d0
!! Dump array
integer, allocatable :: n_kl_dump(:,:)
double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:)
allocate(n_kl_dump(pseudo_kmax,0:pseudo_lmax), v_kl_dump(pseudo_kmax,0:pseudo_lmax), dz_kl_dump(pseudo_kmax,0:pseudo_lmax))
! _
! / _. | _ |
! \_ (_| | (_ |_| |
!
print*, 'Providing the nuclear electron pseudo integrals '
call wall_time(wall_1)
call cpu_time(cpu_1)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
!$OMP num_A,num_B,Z,c,n_pt_in, &
!$OMP n_kl_dump, v_kl_dump, dz_kl_dump, &
!$OMP wall_0,wall_2,thread_num) &
!$OMP SHARED (aux_basis_num_sqrt,aux_basis_prim_num,aux_basis_expo_transp,aux_basis_power,aux_basis_nucl,nucl_coord,aux_basis_coef_transp, &
!$OMP aux_pseudo_integral_non_local,nucl_num,nucl_charge, &
!$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl, &
!$OMP wall_1)
!$OMP DO SCHEDULE (guided)
do j = 1, aux_basis_num_sqrt
num_A = aux_basis_nucl(j)
power_A(1:3)= aux_basis_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, aux_basis_num_sqrt
num_B = aux_basis_nucl(i)
power_B(1:3)= aux_basis_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l=1,aux_basis_prim_num(j)
alpha = aux_basis_expo_transp(l,j)
do m=1,aux_basis_prim_num(i)
beta = aux_basis_expo_transp(m,i)
double precision :: c
c = 0.d0
do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3)
n_kl_dump = pseudo_n_kl(k,1:pseudo_kmax,0:pseudo_lmax)
v_kl_dump = pseudo_v_kl(k,1:pseudo_kmax,0:pseudo_lmax)
dz_kl_dump = pseudo_dz_kl(k,1:pseudo_kmax,0:pseudo_lmax)
c = c + Vpseudo(pseudo_lmax,pseudo_kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center)
enddo
aux_pseudo_integral_non_local(i,j) = aux_pseudo_integral_non_local(i,j) + &
aux_basis_coef_transp(l,j)*aux_basis_coef_transp(m,i)*c
enddo
enddo
enddo
call wall_time(wall_2)
if (thread_num == 0) then
if (wall_2 - wall_0 > 1.d0) then
wall_0 = wall_2
print*, 100.*float(j)/float(aux_basis_num_sqrt), '% in ', &
wall_2-wall_1, 's'
endif
endif
enddo
!$OMP END DO
!$OMP END PARALLEL
deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num,pseudo_grid_size) ]
implicit none
@ -240,9 +17,11 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l
integer :: n_a(3)
double precision :: a(3), c(3), g_a
integer :: i,j,k,l,m,n,p
double precision :: r(pseudo_grid_size), dr, Ulc
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
@ -269,6 +48,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num,-pseudo_lmax:pseudo_l
enddo
enddo
enddo
deallocate(r)
END_PROVIDER

163
src/QmcChem/pseudo.irp.f
View File

@ -3,169 +3,6 @@ subroutine write_pseudopotential
BEGIN_DOC
! Write the pseudo_potential into the EZFIO file
END_DOC
! call ezfio_set_pseudo_pseudo_matrix(pseudo_matrix)
! call ezfio_set_pseudo_ao_pseudo_grid(ao_pseudo_grid)
call ezfio_set_pseudo_mo_pseudo_grid(mo_pseudo_grid)
end
BEGIN_PROVIDER [ double precision, pseudo_matrix, (aux_basis_num_sqrt,aux_basis_num_sqrt) ]
implicit none
BEGIN_DOC
! Pseudo-potential expressed in the basis of ao products
END_DOC
integer :: i,j,k,l
integer :: info, m,n, lwork, lda, ldu, ldvt
integer, allocatable :: iwork(:)
character :: jobz
double precision, allocatable :: a(:,:),work(:)
double precision,allocatable :: U(:,:)
double precision,allocatable :: Vt(:,:)
double precision,allocatable :: S(:), B(:)
jobz = 'A'
m = aux_basis_num
n = aux_basis_num
lda = size(aux_basis_overlap_matrix,1)
ldu = lda
ldvt = lda
lwork = -1
! allocate (A(lda,n), U(ldu,n), Vt(ldvt,n), S(n), work(1), b(n), iwork(8*n))
allocate (A(lda,n), U(ldu,n), Vt(ldvt,n), S(n), work(1), b(n),iwork(1))
work(1) = 1
do i=1,n
do j=1,n
A(i,j) = aux_basis_overlap_matrix(i,j)
enddo
enddo
! call dgesdd(jobz, m, n, A, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
call dgesvd(jobz, jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
lwork = int(work(1))
deallocate(work)
print *, 'Fitting pseudo-potentials'
allocate(work(lwork))
! call dgesdd(jobz, m, n, A, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
call dgesvd(jobz, jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
deallocate(work)
do i=1,n
print *, i, s(i)
enddo
do k=1,n
if (s(k) < 1.d-1) then
s(k) = 0.d0
else
s(k) = 1.d0/s(k)
endif
do m=1,n
Vt(m,k) = S(m) * Vt(m,k)
enddo
enddo
call dgemm('N','N',n,n,n,1.d0,U,lda,Vt,ldvt,0.d0,A,lda)
! do k=1,n
! do l=1,n
! A(k,l) = 0.d0
! do m=1,n
! A(k,l) = A(k,l) + U(k,m) * Vt(m,l)
! enddo
! enddo
do k=1,n
i = aux_basis_idx(1,k)
j = aux_basis_idx(2,k)
b(k) = aux_pseudo_integral(i,j)
enddo
do k=1,n
S(k) = 0.d0
enddo
do l=1,n
do k=1,n
S(k) = S(k) + A(k,l) * b(l)
enddo
enddo
do k=1,aux_basis_num
i = aux_basis_idx(1,k)
j = aux_basis_idx(2,k)
pseudo_matrix(i,j) = S(k)
pseudo_matrix(j,i) = S(k)
enddo
deallocate(a,b,s,iwork,u,vt)
print *, 'Done'
if (info /= 0) then
print *, info
stop 'pseudo fit failed'
endif
END_PROVIDER
!BEGIN_PROVIDER [ double precision, pseudo_matrix, (ao_num,ao_num) ]
! implicit none
! BEGIN_DOC
! ! Pseudo-potential expressed in the basis of ao products
! END_DOC
!
! integer :: i,j,k
! integer :: info, n, lwork, lda, ldb, nrhs
! character :: uplo
! integer, allocatable :: ipiv(:)
! double precision, allocatable :: a(:,:),work(:), b(:)
!
! uplo = 'L'
! n = aux_basis_num
! nrhs = 1
! lda = size(aux_basis_overlap_matrix,1)
! ldb = n
! lwork = -1
!
! print *, 'Fitting pseudo-potentials'
! allocate(work(1),a(lda,n),ipiv(n),b(n))
! work(1) = 1
! do i=1,n
! do j=1,n
! a(i,j) = aux_basis_overlap_matrix(i,j)
! enddo
! enddo
!
! do k=1,n
! i = aux_basis_idx(1,k)
! j = aux_basis_idx(2,k)
! b(k) = ao_pseudo_integral(i,j)
! enddo
! call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info )
! lwork = int(work(1))
! deallocate(work)
!
! allocate(work(lwork))
! call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info )
! deallocate(work,ipiv)
! do k=1,aux_basis_num
! i = aux_basis_idx(1,k)
! j = aux_basis_idx(2,k)
! pseudo_matrix(i,j) = b(k)
! pseudo_matrix(j,i) = b(k)
! enddo
! deallocate(a,b)
!
!print *, 'Done'
! if (info /= 0) then
! print *, info
! stop 'pseudo fit failed'
! endif
!END_PROVIDER