From c8a5cf37cdd26e0d9e6d134ae7f3002e68bb32f2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 23 Sep 2016 16:16:48 +0200 Subject: [PATCH] Working on davidson --- configure | 1 - plugins/Full_CI_ZMQ/fci_zmq.irp.f | 50 ++++++++++++++++-------------- src/Davidson/diagonalization.irp.f | 23 ++++++-------- src/Davidson/u0Hu0.irp.f | 31 ++++++++++++------ 4 files changed, 58 insertions(+), 47 deletions(-) diff --git a/configure b/configure index 2e0d6c49..19016136 100755 --- a/configure +++ b/configure @@ -543,7 +543,6 @@ def recommendation(): print "" print "Finally :" print " ninja" - print " make -C ocaml" print "" print "You can install more plugin with the qp_module.py install command" print "PS : For more info on compiling the code, read the README.md" diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index 0a50a623..dbf8414b 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/plugins/Full_CI_ZMQ/fci_zmq.irp.f @@ -3,12 +3,10 @@ program fci_zmq integer :: i,j,k logical, external :: detEq - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) - integer :: N_st, degree - integer(bit_kind) :: chk + double precision, allocatable :: pt2(:) + integer :: degree - N_st = N_states - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + allocate (pt2(N_states)) pt2 = 1.d0 diag_algorithm = "Lapack" @@ -24,20 +22,22 @@ program fci_zmq call save_wavefunction print *, 'N_det = ', N_det print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 - print *, '-----' + do k=1,N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E+PT2 = ', CI_energy(k) + pt2(k) + print *, '-----' + enddo endif - double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states) double precision :: E_CI_before(N_states) integer :: n_det_before print*,'Beginning the selection ...' - E_CI_before = CI_energy + E_CI_before(1:N_states) = CI_energy(1:N_states) - do while ( (N_det < N_det_max) .and. (maxval(abs(pt2(1:N_st))) > pt2_max) ) + do while ( (N_det < N_det_max) .and. (maxval(abs(pt2(1:N_states))) > pt2_max) ) n_det_before = N_det call ZMQ_selection(max(1024-N_det, N_det), pt2) @@ -56,14 +56,13 @@ program fci_zmq print *, 'N_det = ', N_det print *, 'N_states = ', N_states - do k = 1, N_states - print*,'State ',k - print *, 'PT2 = ', pt2(k) - print *, 'E = ', CI_energy(k) - print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) + do k=1, N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) enddo print *, '-----' - E_CI_before = CI_energy if(N_states.gt.1)then print*,'Variational Energy difference' do i = 2, N_states @@ -76,7 +75,7 @@ program fci_zmq print*,'Delta E = ',E_CI_before(i)+ pt2(i) - (E_CI_before(1) + pt2(1)) enddo endif - E_CI_before = CI_energy + E_CI_before(1:N_states) = CI_energy(1:N_states) call ezfio_set_full_ci_energy(CI_energy) enddo N_det = min(N_det_max,N_det) @@ -86,15 +85,18 @@ program fci_zmq print*,'Last iteration only to compute the PT2' threshold_selectors = 1.d0 threshold_generators = 0.9999d0 - E_CI_before = CI_energy + E_CI_before(1:N_states) = CI_energy(1:N_states) call ZMQ_selection(1, pt2) print *, 'Final step' print *, 'N_det = ', N_det print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', E_CI_before - print *, 'E+PT2 = ', E_CI_before+pt2 - print *, '-----' + do k=1,N_states + print *, 'State', k + print *, 'PT2 = ', pt2 + print *, 'E = ', E_CI_before + print *, 'E+PT2 = ', E_CI_before+pt2 + print *, '-----' + enddo call ezfio_set_full_ci_energy_pt2(E_CI_before+pt2) endif call save_wavefunction diff --git a/src/Davidson/diagonalization.irp.f b/src/Davidson/diagonalization.irp.f index 21b13113..22a77cae 100644 --- a/src/Davidson/diagonalization.irp.f +++ b/src/Davidson/diagonalization.irp.f @@ -505,24 +505,21 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia ! Gram-Schmidt ! ------------ - double precision :: c + double precision :: c(N_st_diag) do k=1,N_st_diag do iter2=1,iter - do l=1,N_st_diag - c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze) - do i=1,sze - U(i,k,iter+1) = U(i,k,iter+1) - c * U(i,l,iter2) - enddo - enddo - enddo - do l=1,k-1 - c = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze) - do i=1,sze - U(i,k,iter+1) = U(i,k,iter+1) - c * U(i,l,iter+1) - enddo + call dgemv('T',sze,N_st_diag,1.d0,U(1,1,iter2),size(U,1), & + U(1,k,iter+1),1,0.d0,c,1) + call dgemv('N',sze,N_st_diag,-1.d0,U(1,1,iter2),size(U,1), & + c,1,1.d0,U(1,k,iter+1),1) enddo + call dgemv('T',sze,N_st_diag,1.d0,U(1,1,iter+1),size(U,1), & + U(1,k,iter+1),1,0.d0,c,1) + call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1), & + c,1,1.d0,U(1,k,iter+1),1) call normalize( U(1,k,iter+1), sze ) enddo + !DEBUG : CHECK OVERLAP !print *, '===' diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index b875ff33..0339c31f 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -43,7 +43,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) double precision, intent(in) :: H_jj(n) integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) double precision :: hij - double precision, allocatable :: vt(:,:) + double precision, allocatable :: vt(:,:), ut(:,:) integer :: i,j,k,l, jj,ii integer :: i0, j0 @@ -52,7 +52,12 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) integer(bit_kind) :: sorted_i(Nint) integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate + integer :: N_st_8 + integer, external :: align_double + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut + + N_st_8 = align_double(N_st) ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -60,15 +65,23 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) PROVIDE ref_bitmask_energy allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) + allocate(ut(N_st_8,n)) + v_0 = 0.d0 + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(i,istate) + enddo + enddo + call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,sze_8) - allocate(vt(sze_8,N_st)) + !$OMP SHARED(n,H_jj,keys_tmp,ut,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) + allocate(vt(N_st_8,n)) Vt = 0.d0 !$OMP DO SCHEDULE(dynamic) @@ -102,8 +115,8 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) if(ext <= 4) then call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) do istate=1,N_st - vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) - vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) + vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) enddo endif enddo @@ -125,8 +138,8 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) if(ext == 4) then call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) do istate=1,N_st - vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) - vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) + vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) enddo end if end do @@ -137,7 +150,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) !$OMP CRITICAL do istate=1,N_st do i=n,1,-1 - v_0(i,istate) = v_0(i,istate) + vt(i,istate) + v_0(i,istate) = v_0(i,istate) + vt(istate,i) enddo enddo !$OMP END CRITICAL @@ -150,7 +163,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) v_0(i,istate) += H_jj(i) * u_0(i,istate) enddo enddo - deallocate (shortcut, sort_idx, sorted, version) + deallocate (shortcut, sort_idx, sorted, version, ut) end BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]