mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-09-01 08:43:46 +02:00
Compare commits
6 Commits
ec1a0e8af4
...
edc592e012
Author | SHA1 | Date | |
---|---|---|---|
edc592e012 | |||
54aaaad227 | |||
2d20035e02 | |||
d8b80f4b55 | |||
cf7e20ff0b | |||
c77ac2598f |
@ -4,3 +4,9 @@ doc: Dressing matrix obtained from DMC
|
|||||||
size: (determinants.n_det)
|
size: (determinants.n_det)
|
||||||
interface: ezfio, provider
|
interface: ezfio, provider
|
||||||
|
|
||||||
|
[dmc_delta_htc]
|
||||||
|
type: double precision
|
||||||
|
doc: Dressing matrix obtained from H_TC
|
||||||
|
size: (determinants.n_det)
|
||||||
|
interface: ezfio
|
||||||
|
|
||||||
|
@ -12,21 +12,26 @@ sys.path.insert(0,QP_PATH)
|
|||||||
from ezfio import ezfio
|
from ezfio import ezfio
|
||||||
|
|
||||||
def read_hamiltonian(inp):
|
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
|
inside = None
|
||||||
h = []
|
h = []
|
||||||
s = []
|
s = []
|
||||||
norm = None
|
norm = None
|
||||||
for line in text.splitlines():
|
for line in text.splitlines():
|
||||||
line = str(line)
|
line = line.decode("utf-8")
|
||||||
print (line)
|
if line == "":
|
||||||
|
continue
|
||||||
if "Psi_norm :" in line:
|
if "Psi_norm :" in line:
|
||||||
norm = float(line.split()[3])
|
norm = float(line.split()[2])
|
||||||
if "]" in line:
|
if "]" in line:
|
||||||
inside = None
|
inside = None
|
||||||
if inside == "H":
|
if inside == "H":
|
||||||
data = line.split()
|
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:
|
elif "Ci_dress" in line:
|
||||||
inside = "H"
|
inside = "H"
|
||||||
h = np.array(h)/norm
|
h = np.array(h)/norm
|
||||||
@ -38,3 +43,4 @@ h = read_hamiltonian(sys.argv[1])
|
|||||||
ezfio.set_file(sys.argv[1])
|
ezfio.set_file(sys.argv[1])
|
||||||
ezfio.set_dmc_dress_dmc_delta_h(h)
|
ezfio.set_dmc_dress_dmc_delta_h(h)
|
||||||
print(h)
|
print(h)
|
||||||
|
|
||||||
|
@ -11,7 +11,7 @@
|
|||||||
double precision, allocatable :: delta(:)
|
double precision, allocatable :: delta(:)
|
||||||
|
|
||||||
allocate(delta(N_det))
|
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)
|
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
|
||||||
|
@ -24,7 +24,8 @@ subroutine run
|
|||||||
double precision, allocatable :: A(:,:)
|
double precision, allocatable :: A(:,:)
|
||||||
double precision, allocatable :: V(:)
|
double precision, allocatable :: V(:)
|
||||||
integer , allocatable :: Vi(:,:)
|
integer , allocatable :: Vi(:,:)
|
||||||
double precision, allocatable :: s
|
double precision :: s
|
||||||
|
PROVIDE ao_num
|
||||||
|
|
||||||
f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
|
f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
|
||||||
if (f == 0_8) then
|
if (f == 0_8) then
|
||||||
@ -37,8 +38,9 @@ subroutine run
|
|||||||
|
|
||||||
if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
|
if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
|
||||||
rc = trexio_read_nucleus_repulsion(f, s)
|
rc = trexio_read_nucleus_repulsion(f, s)
|
||||||
|
call trexio_assert(rc, TREXIO_SUCCESS)
|
||||||
if (rc /= TREXIO_SUCCESS) then
|
if (rc /= TREXIO_SUCCESS) then
|
||||||
print *, irp_here
|
print *, irp_here, rc
|
||||||
print *, 'Error reading nuclear repulsion'
|
print *, 'Error reading nuclear repulsion'
|
||||||
stop -1
|
stop -1
|
||||||
endif
|
endif
|
||||||
@ -100,20 +102,21 @@ subroutine run
|
|||||||
|
|
||||||
! AO 2e integrals
|
! AO 2e integrals
|
||||||
! ---------------
|
! ---------------
|
||||||
|
PROVIDE ao_integrals_map
|
||||||
|
|
||||||
allocate(buffer_i(ao_num**3), buffer_values(ao_num**3))
|
integer*4 :: BUFSIZE
|
||||||
allocate(Vi(4,ao_num**3), V(ao_num**3))
|
BUFSIZE=ao_num**2
|
||||||
|
allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
|
||||||
|
allocate(Vi(4,BUFSIZE), V(BUFSIZE))
|
||||||
|
|
||||||
integer*8 :: offset, icount
|
integer*8 :: offset, icount
|
||||||
|
|
||||||
|
open(unit=104,file='tmp')
|
||||||
offset = 0_8
|
offset = 0_8
|
||||||
icount = 0_8
|
icount = BUFSIZE
|
||||||
rc = TREXIO_SUCCESS
|
rc = TREXIO_SUCCESS
|
||||||
do while (icount == size(V))
|
do while (icount == size(V))
|
||||||
rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
|
rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
|
||||||
if (rc /= TREXIO_SUCCESS) then
|
|
||||||
exit
|
|
||||||
endif
|
|
||||||
do m=1,icount
|
do m=1,icount
|
||||||
i = Vi(1,m)
|
i = Vi(1,m)
|
||||||
j = Vi(2,m)
|
j = Vi(2,m)
|
||||||
@ -122,11 +125,16 @@ subroutine run
|
|||||||
integral = V(m)
|
integral = V(m)
|
||||||
call two_e_integrals_index(i, j, k, l, buffer_i(m) )
|
call two_e_integrals_index(i, j, k, l, buffer_i(m) )
|
||||||
buffer_values(m) = integral
|
buffer_values(m) = integral
|
||||||
|
write(104,'(4(I5,X),2D22.15)') i,j,k,l, integral
|
||||||
enddo
|
enddo
|
||||||
call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
|
call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
|
||||||
offset = offset + icount
|
offset = offset + icount
|
||||||
|
if (rc /= TREXIO_SUCCESS) then
|
||||||
|
exit
|
||||||
|
endif
|
||||||
end do
|
end do
|
||||||
n_integrals = offset
|
n_integrals = offset
|
||||||
|
close(104)
|
||||||
|
|
||||||
call map_sort(ao_integrals_map)
|
call map_sort(ao_integrals_map)
|
||||||
call map_unique(ao_integrals_map)
|
call map_unique(ao_integrals_map)
|
||||||
|
@ -10,4 +10,14 @@ fi
|
|||||||
|
|
||||||
pkg-config --libs trexio > LIB
|
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
|
ln --symbolic ${PWD}/qp_import_trexio.py $QP_ROOT/scripts
|
||||||
|
@ -68,21 +68,13 @@ def write_ezfio(trexio_filename, filename):
|
|||||||
trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_HDF5)
|
trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_HDF5)
|
||||||
|
|
||||||
basis_type = trexio.read_basis_type(trexio_file)
|
basis_type = trexio.read_basis_type(trexio_file)
|
||||||
if basis_type.lower() != "gaussian":
|
# if basis_type.lower() != "gaussian":
|
||||||
raise TypeError
|
# raise TypeError
|
||||||
|
|
||||||
ezfio.set_file(filename)
|
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=' ')
|
print("Nuclei\t\t...\t", end=' ')
|
||||||
|
|
||||||
@ -106,6 +98,23 @@ def write_ezfio(trexio_filename, filename):
|
|||||||
print("OK")
|
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=' ')
|
print("Basis\t\t...\t", end=' ')
|
||||||
|
|
||||||
ezfio.set_basis_basis("Read from TREXIO")
|
ezfio.set_basis_basis("Read from TREXIO")
|
||||||
@ -229,6 +238,7 @@ def write_ezfio(trexio_filename, filename):
|
|||||||
coef.append(coefficient[j])
|
coef.append(coefficient[j])
|
||||||
expo.append(exponent[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_coef(coef)
|
||||||
ezfio.set_ao_basis_ao_expo(expo)
|
ezfio.set_ao_basis_ao_expo(expo)
|
||||||
ezfio.set_ao_basis_ao_basis("Read from TREXIO")
|
ezfio.set_ao_basis_ao_basis("Read from TREXIO")
|
||||||
@ -265,17 +275,21 @@ def write_ezfio(trexio_filename, filename):
|
|||||||
except trexio.Error:
|
except trexio.Error:
|
||||||
pass
|
pass
|
||||||
|
|
||||||
MoMatrix = trexio.read_mo_coefficient(trexio_file)
|
try:
|
||||||
mo_num = trexio.read_mo_num(trexio_file)
|
mo_num = trexio.read_mo_num(trexio_file)
|
||||||
|
ezfio.set_mo_basis_mo_num(mo_num)
|
||||||
|
|
||||||
ezfio.set_mo_basis_mo_num(mo_num)
|
MoMatrix = trexio.read_mo_coefficient(trexio_file)
|
||||||
ezfio.set_mo_basis_mo_coef(MoMatrix)
|
ezfio.set_mo_basis_mo_coef(MoMatrix)
|
||||||
mo_occ = [ 0. for i in range(mo_num) ]
|
|
||||||
for i in range(num_alpha):
|
mo_occ = [ 0. for i in range(mo_num) ]
|
||||||
mo_occ[i] += 1.
|
for i in range(num_alpha):
|
||||||
for i in range(num_beta):
|
mo_occ[i] += 1.
|
||||||
mo_occ[i] += 1.
|
for i in range(num_beta):
|
||||||
ezfio.set_mo_basis_mo_occ(mo_occ)
|
mo_occ[i] += 1.
|
||||||
|
ezfio.set_mo_basis_mo_occ(mo_occ)
|
||||||
|
except:
|
||||||
|
pass
|
||||||
|
|
||||||
print("OK")
|
print("OK")
|
||||||
|
|
||||||
|
@ -8,7 +8,7 @@ then
|
|||||||
exit -1
|
exit -1
|
||||||
fi
|
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
|
# Destroy ONLY the symbolic link for the scripts to be used in the
|
||||||
# ${QP_ROOT}/scripts/ directory.
|
# ${QP_ROOT}/scripts/ directory.
|
||||||
|
17
fnmf/EZFIO.cfg
Normal file
17
fnmf/EZFIO.cfg
Normal 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
|
||||||
|
|
4
fnmf/README.rst
Normal file
4
fnmf/README.rst
Normal file
@ -0,0 +1,4 @@
|
|||||||
|
====
|
||||||
|
fnmf
|
||||||
|
====
|
||||||
|
|
70
fnmf/dmc_data.irp.f
Normal file
70
fnmf/dmc_data.irp.f
Normal 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
76
fnmf/fnmf.irp.f
Normal 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
74
fnmf/non_hermit.irp.f
Normal 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
|
||||||
|
|
Loading…
Reference in New Issue
Block a user