From 9548cf14d4d544bf2f3dc5a8d9ece8688f75da01 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 21 Jun 2015 23:08:18 +0200 Subject: [PATCH] Work on QMC=Chem interface --- config/gfortran.cfg | 2 +- scripts/compilation/qp_create_ninja.py | 2 +- src/QmcChem/pot_ao_pseudo_ints.irp.f | 228 +------------------------ src/QmcChem/pseudo.irp.f | 163 ------------------ 4 files changed, 6 insertions(+), 389 deletions(-) diff --git a/config/gfortran.cfg b/config/gfortran.cfg index d84d8486..12cd6e0c 100644 --- a/config/gfortran.cfg +++ b/config/gfortran.cfg @@ -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 diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index 9277b276..91f42c81 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -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)), diff --git a/src/QmcChem/pot_ao_pseudo_ints.irp.f b/src/QmcChem/pot_ao_pseudo_ints.irp.f index 0aa2a98c..3b655b00 100644 --- a/src/QmcChem/pot_ao_pseudo_ints.irp.f +++ b/src/QmcChem/pot_ao_pseudo_ints.irp.f @@ -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 diff --git a/src/QmcChem/pseudo.irp.f b/src/QmcChem/pseudo.irp.f index abf4a5b8..13981002 100644 --- a/src/QmcChem/pseudo.irp.f +++ b/src/QmcChem/pseudo.irp.f @@ -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 -