scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-12-22 12:23:37 +01:00
Code Issues Releases Wiki Activity

Merge branch 'master' of gitlab.com:scemama/qp_plugins_scemama

Browse Source
This commit is contained in:
Anthony Scemama 2023-03-09 10:52:12 +01:00
commit 533026d3cb
11 changed files with 299 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
devel/cc/denominators.irp.f
View File

@ -21,6 +21,7 @@ BEGIN_PROVIDER [ double precision, delta_OV, (spin_occ_num, spin_vir_num) ]
do i=1,spin_occ_num
delta_OV(i,a) = eO(i) - eV(a)
! delta_OV(i,a) = min(eO(i) - eV(a), -1.d-3)
if (delta_OV(i,a) == 0.d0) delta_OV(i,a) = 1.d-1
enddo
enddo
@ -39,6 +40,7 @@ BEGIN_PROVIDER [ double precision, delta_OOVV, (spin_occ_num,spin_occ_num,spin_v
do j=1,spin_occ_num
do i=1,spin_occ_num
delta_OOVV(i,j,a,b) = delta_OV(i,a) + delta_OV(j,b)
if (delta_OOVV(i,j,a,b) == 0.d0) delta_OOVV(i,j,a,b) = 1.d-1
enddo
enddo
enddo
@ -60,6 +62,7 @@ BEGIN_PROVIDER [ double precision, delta_OOOVVV, (spin_occ_num,spin_occ_num,spin
do j=1,spin_occ_num
do i=1,spin_occ_num
delta_OOOVVV(i,j,k,a,b,c) = delta_OV(i,a) + delta_OV(j,b) + delta_OV(k,c)
if (delta_OOOVVV(i,j,k,a,b,c) == 0.d0) delta_OOOVVV(i,j,k,a,b,c) = 1.d-1
enddo
enddo
enddo

6
devel/dmc_dress/EZFIO.cfg
View File

@ -4,3 +4,9 @@ doc: Dressing matrix obtained from DMC
size: (determinants.n_det)
interface: ezfio, provider
[dmc_delta_htc]
type: double precision
doc: Dressing matrix obtained from H_TC
size: (determinants.n_det)
interface: ezfio

16
devel/dmc_dress/dress.py
View File

@ -12,21 +12,26 @@ sys.path.insert(0,QP_PATH)
from ezfio import ezfio
def read_hamiltonian(inp):
text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
# text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
text = subprocess.run(["cat", "OUT"], capture_output=True).stdout
inside = None
h = []
s = []
norm = None
for line in text.splitlines():
line = str(line)
print (line)
line = line.decode("utf-8")
if line == "":
continue
if "Psi_norm :" in line:
norm = float(line.split()[3])
norm = float(line.split()[2])
if "]" in line:
inside = None
if inside == "H":
data = line.split()
h.append(float(data[3]))
# h.append(float(data[2])) #ave
h.append(max(float(data[2]) - float(data[4]),0.)) #diff
# if float(data[2]) - float(data[4])>0.: h.append(float(data[2]))
# else: h.append(0.)
elif "Ci_dress" in line:
inside = "H"
h = np.array(h)/norm
@ -38,3 +43,4 @@ h = read_hamiltonian(sys.argv[1])
ezfio.set_file(sys.argv[1])
ezfio.set_dmc_dress_dmc_delta_h(h)
print(h)

37
devel/dmc_dress/dressing_vector.irp.f
View File

@ -11,7 +11,7 @@
double precision, allocatable :: delta(:)
allocate(delta(N_det))
delta(1:N_det) = dmc_delta_h(1:N_det)
delta (1:N_det) = dmc_delta_h(1:N_det) + dmc_delta_htc(1:N_det)
call dset_order(delta,psi_bilinear_matrix_order_reverse,N_det)
@ -31,3 +31,38 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, dmc_delta_htc , (n_det) ]
implicit none
BEGIN_DOC
! Dressing matrix obtained from H_TC
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
if (size(dmc_delta_htc) == 0) return
call ezfio_has_dmc_dress_dmc_delta_htc(has)
if (has) then
write(6,'(A)') '.. >>>>> [ IO READ: dmc_delta_htc ] <<<<< ..'
call ezfio_get_dmc_dress_dmc_delta_htc(dmc_delta_htc)
else
dmc_delta_htc(:) = 0.d0
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( dmc_delta_htc, (n_det), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read dmc_delta_htc with MPI'
endif
IRP_ENDIF
call write_time(6)
END_PROVIDER

3
devel/trexio/import_trexio_integrals.irp.f
View File

@ -111,7 +111,6 @@ subroutine run
integer*8 :: offset, icount
! open(unit=104,file='tmp')
offset = 0_8
icount = BUFSIZE
rc = TREXIO_SUCCESS
@ -125,7 +124,6 @@ subroutine run
integral = V(m)
call two_e_integrals_index(i, j, k, l, buffer_i(m) )
buffer_values(m) = integral
! write(104,'(4(I5,X),2D22.15)') i,j,k,l, integral
enddo
call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
offset = offset + icount
@ -134,7 +132,6 @@ subroutine run
endif
end do
n_integrals = offset
close(104)
call map_sort(ao_integrals_map)
call map_unique(ao_integrals_map)

17
fnmf/EZFIO.cfg Normal file
View File

@ -0,0 +1,17 @@
[E_dmc]
type: double precision
doc: DMC energy
interface: ezfio, provider
[h_dmc]
type: double precision
doc: Dressing matrix obtained from DMC
size: (determinants.n_det)
interface: ezfio, provider
[s_dmc]
type: double precision
doc: Dressing matrix obtained from H_TC
size: (determinants.n_det)
interface: ezfio, provider

2
fnmf/NEED Normal file
View File

@ -0,0 +1,2 @@
determinants
davidson_undressed

4
fnmf/README.rst Normal file
View File

@ -0,0 +1,4 @@
====
fnmf
====

70
fnmf/dmc_data.irp.f Normal file
View File

@ -0,0 +1,70 @@
BEGIN_PROVIDER [ double precision, h_dmc_row , (N_det) ]
&BEGIN_PROVIDER [ double precision, s_dmc_row , (N_det) ]
implicit none
BEGIN_DOC
! Data sampled with QMC=Chem
END_DOC
h_dmc_row(:) = h_dmc(:)
s_dmc_row(:) = s_dmc(:)
call dset_order(h_dmc_row,psi_bilinear_matrix_order_reverse,N_det)
call dset_order(s_dmc_row,psi_bilinear_matrix_order_reverse,N_det)
! integer :: i
! do i=1,N_det
! s_dmc_row(i) = psi_coef(i,1)
! call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
! N_det, 1, h_dmc_row(i) )
! enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, mat_size ]
implicit none
BEGIN_DOC
! Size of the matrices
END_DOC
mat_size = N_det+1
END_PROVIDER
BEGIN_PROVIDER [ double precision, H_dmc_mat, (mat_size, mat_size) ]
implicit none
BEGIN_DOC
! Hamiltonian extended with DMC data
END_DOC
integer :: i,j
do j=1,N_det
do i=1,N_det
call i_h_j(psi_det(1,1,i), psi_det(1,1,j), N_int, H_dmc_mat(i,j))
enddo
enddo
do i=1,N_det
call i_h_psi(psi_det(1,1,i), psi_det, psi_coef, N_int, N_det, &
N_det, 1, H_dmc_mat(i,N_det+1) )
H_dmc_mat(N_det+1,i) = h_dmc_row(i)
enddo
H_dmc_mat(mat_size,mat_size) = E_dmc - nuclear_repulsion
END_PROVIDER
BEGIN_PROVIDER [ double precision, S_dmc_mat, (mat_size, mat_size) ]
implicit none
BEGIN_DOC
! Overlap matrix extended with DMC data
END_DOC
integer :: i,j
S_dmc_mat = 0.d0
do i=1,mat_size
S_dmc_mat(i,i) = 1.d0
enddo
do i=1,N_det
S_dmc_mat(i,N_det+1) = psi_coef(i,1)
S_dmc_mat(N_det+1,i) = S_dmc_row(i)
enddo
END_PROVIDER

76
fnmf/fnmf.irp.f Normal file
View File

@ -0,0 +1,76 @@
program fnmf
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
read_wf = .True.
TOUCH read_wf
call run
end
subroutine run
implicit none
integer :: i, n_real
double precision, allocatable :: beta(:), vr(:,:), vl(:,:)
double precision, allocatable :: htmp(:,:), stmp(:,:)
allocate(Htmp(N_det,N_det), Stmp(N_det,N_det))
htmp(:, :) = H_dmc_mat(1:N_det, 1:N_det)
stmp(:, :) = S_dmc_mat(1:N_det, 1:N_det)
htmp(1, 1) = H_dmc_mat(N_det+1, N_det+1)
stmp(1, 1) = S_dmc_mat(N_det+1, N_det+1)
htmp(2:N_det, 1) = H_dmc_mat(2:N_det, N_det+1)
stmp(2:N_det, 1) = S_dmc_mat(2:N_det, N_det+1)
htmp(1, 2:N_det) = H_dmc_mat(N_det+1, 2:N_det)
stmp(1, 2:N_det) = S_dmc_mat(N_det+1, 2:N_det)
allocate(beta(N_det), vr(N_det,N_det), vl(N_det, N_det))
call lapack_g_non_sym_real(N_det, Htmp, size(Htmp,1), &
Stmp, size(Stmp,1), beta, n_real, vl, size(vl,1), vr, size(vr,1))
vl(:,1) = psi_coef(:,1)
psi_coef(:,1) *= vr(1,1)
do i=2,N_det
psi_coef(i,1) += vr(i,1)
end do
double precision, external :: dnrm2, ddot
double precision :: norm
norm = dnrm2(N_det, psi_coef, 1)
psi_coef(:,1) *= 1.d0/norm
SOFT_TOUCH psi_coef
print *, '-------------------------------------------------'
print *, ' psi_old psi_new'
print *, '-------------------------------------------------'
do i=1,N_det
print '(2(E16.8,X))', vl(i,1), psi_coef(i,1)
enddo
print *, '-------------------------------------------------'
print *, ''
print *, 'Overlap : ', ddot(N_det, psi_coef, 1, vl, 1)
print *, 'DMC energy : ', htmp(1,1) + nuclear_repulsion
print *, 'Updated energy : ', beta(1) + nuclear_repulsion
print *, 'H'
do i=1,N_det
print *, htmp(i, 1:N_det)
enddo
print *, 'S'
do i=1,N_det
print *, stmp(i, 1:N_det)
enddo
call save_wavefunction
end

74
fnmf/non_hermit.irp.f Normal file
View File

@ -0,0 +1,74 @@
subroutine lapack_g_non_sym_real(n, H, LDH, S, LDS, beta, &
n_real, vl, LDVL, vr, LDVR)
implicit none
integer, intent(in) :: n, LDH, LDS, LDVL, LDVR
double precision, intent(in) :: H(LDH,n), S(LDS,n)
double precision, intent(out) :: VL(LDVL,n), VR(LDVR,n), beta(n)
integer, intent(out) :: n_real
integer :: lwork, info, i,j
double precision, allocatable :: work(:), Htmp(:,:), Stmp(:,:)
double precision, allocatable :: alphar(:), alphai(:), vltmp(:,:), vrtmp(:,:)
integer, allocatable :: iorder(:), list_good(:)
lwork = -1
allocate(work(1), alphar(n), alphai(n), vltmp(n,n), vrtmp(n,n), &
Htmp(n,n), Stmp(n,n), list_good(n))
Htmp(1:n,1:n) = H(1:n,1:n)
Stmp(1:n,1:n) = S(1:n,1:n)
call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
call dggev('V', 'V', n, Htmp, size(Htmp,1), Stmp, size(Stmp,1), alphar, alphai, beta, &
vltmp, size(vltmp,1), vrtmp, size(vrtmp,1), work, lwork, info)
deallocate(work)
if (info /= 0) then
stop 'DGGEV Diagonalization failed'
endif
allocate(iorder(n))
n_real = 0
do i=1,n
if (dabs(alphai(i)) < 1.d-10) then
n_real += 1
list_good(n_real) = i
else
alphar(i) = dble(huge(1.0))
endif
enddo
do i=1,n
if (dabs(beta(i)/(alphar(i)+1.d-12)) < 1.d-10) then
beta(i) = 0.d0
else
beta(i) = alphar(i)/beta(i)
endif
enddo
do i=1,n_real
iorder(i) = i
do j=1,n
vrtmp(j,i) = vrtmp(j,list_good(i))
vltmp(j,i) = vltmp(j,list_good(i))
end do
end do
call dsort(beta, iorder, n_real)
do i=1,n_real
do j=1,n
vr(j,i) = vrtmp(j,iorder(i))
vl(j,i) = vltmp(j,iorder(i))
end do
end do
deallocate(vrtmp, vltmp, iorder, Htmp, Stmp)
end