mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-09-01 08:43:46 +02:00
Compare commits
No commits in common. "edc592e012c1ba2da66b0806937413d876978549" and "ec1a0e8af4c2fdebc0883b01aeb2c4b4adeb6405" have entirely different histories.
edc592e012
...
ec1a0e8af4
@ -4,9 +4,3 @@ 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
|
||||
|
||||
|
||||
@ -12,26 +12,21 @@ 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(["cat", "OUT"], capture_output=True).stdout
|
||||
text = subprocess.run(["qmcchem", "result", inp], capture_output=True).stdout
|
||||
inside = None
|
||||
h = []
|
||||
s = []
|
||||
norm = None
|
||||
for line in text.splitlines():
|
||||
line = line.decode("utf-8")
|
||||
if line == "":
|
||||
continue
|
||||
line = str(line)
|
||||
print (line)
|
||||
if "Psi_norm :" in line:
|
||||
norm = float(line.split()[2])
|
||||
norm = float(line.split()[3])
|
||||
if "]" in line:
|
||||
inside = None
|
||||
if inside == "H":
|
||||
data = line.split()
|
||||
# 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.)
|
||||
h.append(float(data[3]))
|
||||
elif "Ci_dress" in line:
|
||||
inside = "H"
|
||||
h = np.array(h)/norm
|
||||
@ -43,4 +38,3 @@ h = read_hamiltonian(sys.argv[1])
|
||||
ezfio.set_file(sys.argv[1])
|
||||
ezfio.set_dmc_dress_dmc_delta_h(h)
|
||||
print(h)
|
||||
|
||||
|
||||
@ -11,7 +11,7 @@
|
||||
double precision, allocatable :: delta(:)
|
||||
|
||||
allocate(delta(N_det))
|
||||
delta (1:N_det) = dmc_delta_h(1:N_det) + dmc_delta_htc(1:N_det)
|
||||
delta(1:N_det) = dmc_delta_h(1:N_det)
|
||||
|
||||
call dset_order(delta,psi_bilinear_matrix_order_reverse,N_det)
|
||||
|
||||
@ -31,38 +31,3 @@ 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
|
||||
|
||||
@ -24,8 +24,7 @@ subroutine run
|
||||
double precision, allocatable :: A(:,:)
|
||||
double precision, allocatable :: V(:)
|
||||
integer , allocatable :: Vi(:,:)
|
||||
double precision :: s
|
||||
PROVIDE ao_num
|
||||
double precision, allocatable :: s
|
||||
|
||||
f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
|
||||
if (f == 0_8) then
|
||||
@ -38,9 +37,8 @@ subroutine run
|
||||
|
||||
if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
|
||||
rc = trexio_read_nucleus_repulsion(f, s)
|
||||
call trexio_assert(rc, TREXIO_SUCCESS)
|
||||
if (rc /= TREXIO_SUCCESS) then
|
||||
print *, irp_here, rc
|
||||
print *, irp_here
|
||||
print *, 'Error reading nuclear repulsion'
|
||||
stop -1
|
||||
endif
|
||||
@ -102,21 +100,20 @@ subroutine run
|
||||
|
||||
! AO 2e integrals
|
||||
! ---------------
|
||||
PROVIDE ao_integrals_map
|
||||
|
||||
integer*4 :: BUFSIZE
|
||||
BUFSIZE=ao_num**2
|
||||
allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
|
||||
allocate(Vi(4,BUFSIZE), V(BUFSIZE))
|
||||
allocate(buffer_i(ao_num**3), buffer_values(ao_num**3))
|
||||
allocate(Vi(4,ao_num**3), V(ao_num**3))
|
||||
|
||||
integer*8 :: offset, icount
|
||||
|
||||
open(unit=104,file='tmp')
|
||||
offset = 0_8
|
||||
icount = BUFSIZE
|
||||
icount = 0_8
|
||||
rc = TREXIO_SUCCESS
|
||||
do while (icount == size(V))
|
||||
rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
|
||||
if (rc /= TREXIO_SUCCESS) then
|
||||
exit
|
||||
endif
|
||||
do m=1,icount
|
||||
i = Vi(1,m)
|
||||
j = Vi(2,m)
|
||||
@ -125,16 +122,11 @@ 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
|
||||
if (rc /= TREXIO_SUCCESS) then
|
||||
exit
|
||||
endif
|
||||
end do
|
||||
n_integrals = offset
|
||||
close(104)
|
||||
|
||||
call map_sort(ao_integrals_map)
|
||||
call map_unique(ao_integrals_map)
|
||||
|
||||
@ -10,14 +10,4 @@ fi
|
||||
|
||||
pkg-config --libs trexio > LIB
|
||||
|
||||
scripts_list="qp_import_trexio.py"
|
||||
|
||||
# Destroy ONLY the symbolic link for the scripts to be used in the
|
||||
# ${QP_ROOT}/scripts/ directory.
|
||||
for i in $scripts_list
|
||||
do
|
||||
find ${QP_ROOT}/scripts/$i -type l -delete
|
||||
done
|
||||
|
||||
# Create symlink in scripts
|
||||
ln --symbolic ${PWD}/qp_import_trexio.py $QP_ROOT/scripts
|
||||
|
||||
@ -68,13 +68,21 @@ def write_ezfio(trexio_filename, filename):
|
||||
trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_HDF5)
|
||||
|
||||
basis_type = trexio.read_basis_type(trexio_file)
|
||||
# if basis_type.lower() != "gaussian":
|
||||
# raise TypeError
|
||||
if basis_type.lower() != "gaussian":
|
||||
raise TypeError
|
||||
|
||||
ezfio.set_file(filename)
|
||||
ezfio.set_trexio_trexio_file(trexio_filename)
|
||||
|
||||
|
||||
print("Electrons\t...\t", end=' ')
|
||||
|
||||
num_alpha = trexio.read_electron_up_num(trexio_file)
|
||||
num_beta = trexio.read_electron_dn_num(trexio_file)
|
||||
ezfio.set_electrons_elec_alpha_num(num_alpha)
|
||||
ezfio.set_electrons_elec_beta_num(num_beta)
|
||||
|
||||
print("OK")
|
||||
|
||||
|
||||
print("Nuclei\t\t...\t", end=' ')
|
||||
|
||||
@ -98,23 +106,6 @@ def write_ezfio(trexio_filename, filename):
|
||||
print("OK")
|
||||
|
||||
|
||||
print("Electrons\t...\t", end=' ')
|
||||
|
||||
try:
|
||||
num_beta = trexio.read_electron_dn_num(trexio_file)
|
||||
except:
|
||||
num_beta = sum(charge)//2
|
||||
|
||||
try:
|
||||
num_alpha = trexio.read_electron_up_num(trexio_file)
|
||||
except:
|
||||
num_alpha = sum(charge) - num_beta
|
||||
|
||||
ezfio.set_electrons_elec_alpha_num(num_alpha)
|
||||
ezfio.set_electrons_elec_beta_num(num_beta)
|
||||
|
||||
print("OK")
|
||||
|
||||
print("Basis\t\t...\t", end=' ')
|
||||
|
||||
ezfio.set_basis_basis("Read from TREXIO")
|
||||
@ -238,7 +229,6 @@ def write_ezfio(trexio_filename, filename):
|
||||
coef.append(coefficient[j])
|
||||
expo.append(exponent[j])
|
||||
|
||||
ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
|
||||
ezfio.set_ao_basis_ao_coef(coef)
|
||||
ezfio.set_ao_basis_ao_expo(expo)
|
||||
ezfio.set_ao_basis_ao_basis("Read from TREXIO")
|
||||
@ -275,21 +265,17 @@ def write_ezfio(trexio_filename, filename):
|
||||
except trexio.Error:
|
||||
pass
|
||||
|
||||
try:
|
||||
mo_num = trexio.read_mo_num(trexio_file)
|
||||
ezfio.set_mo_basis_mo_num(mo_num)
|
||||
MoMatrix = trexio.read_mo_coefficient(trexio_file)
|
||||
mo_num = trexio.read_mo_num(trexio_file)
|
||||
|
||||
MoMatrix = trexio.read_mo_coefficient(trexio_file)
|
||||
ezfio.set_mo_basis_mo_coef(MoMatrix)
|
||||
|
||||
mo_occ = [ 0. for i in range(mo_num) ]
|
||||
for i in range(num_alpha):
|
||||
mo_occ[i] += 1.
|
||||
for i in range(num_beta):
|
||||
mo_occ[i] += 1.
|
||||
ezfio.set_mo_basis_mo_occ(mo_occ)
|
||||
except:
|
||||
pass
|
||||
ezfio.set_mo_basis_mo_num(mo_num)
|
||||
ezfio.set_mo_basis_mo_coef(MoMatrix)
|
||||
mo_occ = [ 0. for i in range(mo_num) ]
|
||||
for i in range(num_alpha):
|
||||
mo_occ[i] += 1.
|
||||
for i in range(num_beta):
|
||||
mo_occ[i] += 1.
|
||||
ezfio.set_mo_basis_mo_occ(mo_occ)
|
||||
|
||||
print("OK")
|
||||
|
||||
|
||||
@ -8,7 +8,7 @@ then
|
||||
exit -1
|
||||
fi
|
||||
|
||||
scripts_list="qp_import_trexio.py"
|
||||
scripts_list="qp_import_trexio.py "
|
||||
|
||||
# Destroy ONLY the symbolic link for the scripts to be used in the
|
||||
# ${QP_ROOT}/scripts/ directory.
|
||||
|
||||
@ -1,17 +0,0 @@
|
||||
[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
|
||||
|
||||
@ -1,4 +0,0 @@
|
||||
====
|
||||
fnmf
|
||||
====
|
||||
|
||||
@ -1,70 +0,0 @@
|
||||
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
|
||||
|
||||
|
||||
|
||||
@ -1,76 +0,0 @@
|
||||
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
|
||||
|
||||
@ -1,74 +0,0 @@
|
||||
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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user