1
0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-12-22 20:34:05 +01:00

Compare commits

..

5 Commits

Author SHA1 Message Date
16d5b14f36 Added svdwf module 2020-09-30 20:44:53 +02:00
04c75ab70c Fixed basis bug 2020-09-24 10:07:25 +02:00
101012938f Removed H matrix elements in champ 2020-08-26 10:53:46 +02:00
164bf79053 Fixed Python3 2020-05-13 00:54:02 +02:00
1568f6cd20 Python3 2020-05-13 00:10:32 +02:00
8 changed files with 379 additions and 160 deletions

View File

@ -96,7 +96,7 @@ subroutine get_fock_matrix_alpha(det,F)
integer :: i,j,k integer :: i,j,k
F(:,:) = fock_op_cshell_ref_bitmask(:,:) F(:,:) = fock_operator_closed_shell_ref_bitmask(:,:)
end end
@ -112,7 +112,7 @@ subroutine get_fock_matrix_beta(det,F)
integer :: i,j,k integer :: i,j,k
F(:,:) = fock_op_cshell_ref_bitmask(:,:) F(:,:) = fock_operator_closed_shell_ref_bitmask(:,:)
end end

59
devel/svdwf/.gitignore vendored Normal file
View File

@ -0,0 +1,59 @@
IRPF90_temp/
IRPF90_man/
build.ninja
irpf90.make
ezfio_interface.irp.f
irpf90_entities
tags
Makefile
ao_basis
ao_one_e_ints
ao_two_e_erf_ints
ao_two_e_ints
aux_quantities
becke_numerical_grid
bitmask
cis
cisd
cipsi
davidson
davidson_dressed
davidson_undressed
density_for_dft
determinants
dft_keywords
dft_utils_in_r
dft_utils_one_e
dft_utils_two_body
dressing
dummy
electrons
ezfio_files
fci
generators_cas
generators_full
hartree_fock
iterations
kohn_sham
kohn_sham_rs
mo_basis
mo_guess
mo_one_e_ints
mo_two_e_erf_ints
mo_two_e_ints
mpi
mrpt_utils
nuclei
perturbation
pseudo
psiref_cas
psiref_utils
scf_utils
selectors_cassd
selectors_full
selectors_utils
single_ref_method
slave
tools
utils
zmq

1
devel/svdwf/NEED Normal file
View File

@ -0,0 +1 @@
determinants

4
devel/svdwf/README.rst Normal file
View File

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

View File

@ -0,0 +1,101 @@
program svdwf
implicit none
BEGIN_DOC
! Make the SVD of the alpha-beta wave function and print singular values.
END_DOC
read_wf = .True.
TOUCH read_wf
call run()
end
subroutine run
implicit none
include 'constants.include.F'
double precision, allocatable :: U(:,:), V(:,:), D(:), A(:,:)
integer :: i, j, k, l, q, r, m, n, iter
double precision,allocatable :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:)
double precision :: r1,r2
m = n_det_alpha_unique
n = n_det_beta_unique
r = min(1000,n)
allocate(Z(m,r))
! Z(m,r) = A(m,n).P(n,r)
Z(:,:) = 0.d0
do l=1,r
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
call random_number(r1)
call random_number(r2)
r1 = dsqrt(-2.d0*dlog(r1))
r2 = dtwo_pi*r2
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1*dcos(r2)
enddo
enddo
! Power iterations
allocate(P(n,r))
do iter=1,20
! P(n,r) = At(n,m).Z(m,r)
P(:,:) = 0.d0
do l=1,r
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
P(j,l) = P(j,l) + psi_bilinear_matrix_values(k,1) * Z(i,l)
enddo
enddo
Z(:,:) = 0.d0
do l=1,r
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * P(j,l)
enddo
enddo
! Compute QR
call ortho_qr(Z,size(Z,1),m,r)
enddo
! Y(r,n) = Zt(r,m).A(m,n)
allocate(Yt(n,r))
Yt(:,:) = 0.d0
do l=1,r
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
Yt(j,l) = Yt(j,l) + Z(i,l) * psi_bilinear_matrix_values(k,1)
enddo
enddo
allocate(D(r),V(n,r), UYt(r,r))
call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r)
deallocate(Yt)
! U(m,r) = Z(m,r).UY(r,r)
allocate(U(m,r))
call dgemm('N','T',m,r,r,1.d0,Z,size(Z,1),UYt,size(UYt,1),0.d0,U,size(U,1))
deallocate(UYt,Z)
do i=1,r
print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
if (D(i) < 1.d-15) then
k = i
exit
endif
enddo
print *, 'threshold: ', 2.858 * D(k/2)
do i=1,m
print '(I6,4(X,F12.8))', i, U(i,1:4)
enddo
print *, ''
do i=1,n
print '(I6,4(X,F12.8))', i, V(i,1:4)
enddo
deallocate(U,D,V)
end

47
devel/svdwf/svdwf.irp.f Normal file
View File

@ -0,0 +1,47 @@
program svdwf
implicit none
BEGIN_DOC
! TODO : Make the SVD of the alpha-beta wave function and print singular values.
END_DOC
read_wf = .True.
TOUCH read_wf
call run()
end
subroutine run
implicit none
double precision, allocatable :: U(:,:), Vt(:,:), D(:), A(:,:)
integer :: i, j, k, p, q
allocate( A (n_det_alpha_unique, n_det_beta_unique), &
U (n_det_alpha_unique, n_det_alpha_unique), &
Vt(n_det_beta_unique, n_det_beta_unique), &
D(max(n_det_beta_unique,n_det_alpha_unique)) )
A = 0.D0
do k=1,N_det
i = psi_bilinear_matrix_rows(k)
j = psi_bilinear_matrix_columns(k)
A(i,j) = psi_bilinear_matrix_values(k,1)
enddo
call randomized_svd(A, size(A,1), &
U, size(U,1), D, Vt, size(Vt,1), n_det_alpha_unique, n_det_beta_unique, &
6,1000)
do i=1,n_det_beta_unique
print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
if (D(i) < 1.d-15) then
k = i
exit
endif
enddo
print *, 'threshold: ', 2.858 * D(k/2)
do i=1,n_det_alpha_unique
print '(I6,4(X,F12.8))', i, U(i,1:4)
enddo
print *, ''
do i=1,n_det_beta_unique
print '(I6,4(X,F12.8))', i, Vt(1:4,i)
enddo
end

View File

@ -2,7 +2,7 @@
# #
# Modified from the QMCPACK interface developed by @tapplencourt and @abenali # Modified from the QMCPACK interface developed by @tapplencourt and @abenali
print "#QP -> CHAMP" print("#QP -> CHAMP")
# ___ # ___
# | ._ o _|_ # | ._ o _|_
@ -19,7 +19,7 @@ ezfio.set_file(ezfio_path)
do_pseudo = ezfio.get_pseudo_do_pseudo() do_pseudo = ezfio.get_pseudo_do_pseudo()
if do_pseudo: if do_pseudo:
print "do_pseudo True" print("do_pseudo True")
from qp_path import QP_ROOT from qp_path import QP_ROOT
l_ele_path = os.path.join(QP_ROOT,"data","list_element.txt") l_ele_path = os.path.join(QP_ROOT,"data","list_element.txt")
@ -28,9 +28,9 @@ if do_pseudo:
l_element_raw = data_raw.split("\n") l_element_raw = data_raw.split("\n")
l_element = [element_raw.split() for element_raw in l_element_raw] l_element = [element_raw.split() for element_raw in l_element_raw]
d_z = dict((abr, z) for (z, abr, ele, _) in filter(lambda x: x != [], l_element) ) d_z = dict((abr, z) for (z, abr, ele, _) in [x for x in l_element if x != []] )
else: else:
print "do_pseudo False" print("do_pseudo False")
try: try:
n_det = ezfio.get_determinants_n_det() n_det = ezfio.get_determinants_n_det()
@ -38,9 +38,9 @@ except IOError:
n_det = 1 n_det = 1
if n_det == 1: if n_det == 1:
print "multi_det False" print("multi_det False")
else: else:
print "multi_det True" print("multi_det True")
# #
# |\/| o _ _ # |\/| o _ _
@ -52,17 +52,17 @@ def list_to_string(l):
ao_num = ezfio.get_ao_basis_ao_num() ao_num = ezfio.get_ao_basis_ao_num()
print "ao_num", ao_num print("ao_num", ao_num)
mo_num = ezfio.get_mo_basis_mo_num() mo_num = ezfio.get_mo_basis_mo_num()
print "mo_num", mo_num print("mo_num", mo_num)
alpha = ezfio.get_electrons_elec_alpha_num() alpha = ezfio.get_electrons_elec_alpha_num()
beta = ezfio.get_electrons_elec_beta_num() beta = ezfio.get_electrons_elec_beta_num()
print "elec_alpha_num", alpha print("elec_alpha_num", alpha)
print "elec_beta_num", beta print("elec_beta_num", beta)
print "elec_tot_num", alpha + beta print("elec_tot_num", alpha + beta)
print "spin_multiplicity", (alpha - beta) + 1 print("spin_multiplicity", (alpha - beta) + 1)
l_label = ezfio.get_nuclei_nucl_label() l_label = ezfio.get_nuclei_nucl_label()
l_charge = ezfio.get_nuclei_nucl_charge() l_charge = ezfio.get_nuclei_nucl_charge()
@ -70,19 +70,19 @@ l_coord = ezfio.get_nuclei_nucl_coord()
l_coord_str = [list_to_string(i) for i in zip(*l_coord)] l_coord_str = [list_to_string(i) for i in zip(*l_coord)]
print "nucl_num", len(l_label) print("nucl_num", len(l_label))
# _ # _
# / _ _ ._ _| # / _ _ ._ _|
# \_ (_) (_) | (_| # \_ (_) (_) | (_|
# #
print "Atomic coord in Bohr" print("Atomic coord in Bohr")
for i, t in enumerate(zip(l_label, l_charge, l_coord_str)): for i, t in enumerate(zip(l_label, l_charge, l_coord_str)):
t_1 = d_z[t[0]] if do_pseudo else t[1] t_1 = d_z[t[0]] if do_pseudo else t[1]
t_new = [t[0],t_1,t[2]] t_new = [t[0],t_1,t[2]]
print list_to_string(t_new) print(list_to_string(t_new))
# #
# Call externet process to get the sysmetry # Call externet process to get the sysmetry
@ -93,14 +93,22 @@ process = subprocess.Popen(
stdout=subprocess.PIPE) stdout=subprocess.PIPE)
out, err = process.communicate() out, err = process.communicate()
basis_raw, sym_raw, _ = out.split("\n\n\n") basis_raw, sym_raw, _ = out.decode().split("\n\n\n")
basis_split = basis_raw.split('\n')
# _ __ # _ __
# |_) _. _ o _ (_ _ _|_ # |_) _. _ o _ (_ _ _|_
# |_) (_| _> | _> __) (/_ |_ # |_) (_| _> | _> __) (/_ |_
# #
basis_without_header = "\n".join(basis_raw.split("\n")[11:]) beginning = 0
for x in basis_split:
if x.startswith("Basis set"):
break
beginning += 1
beginning+=2
basis_without_header = "\n".join(basis_split[beginning:])
import re import re
l_basis_raw = re.split('\n\s*\n', basis_without_header) l_basis_raw = re.split('\n\s*\n', basis_without_header)
@ -118,9 +126,9 @@ for i, (a,b) in enumerate(zip(l_label,l_basis_raw)):
else: else:
continue continue
print "BEGIN_BASIS_SET\n" print("BEGIN_BASIS_SET\n")
print "\n\n".join(l_basis_clean) print("\n\n".join(l_basis_clean))
print "END_BASIS_SET" print("END_BASIS_SET")
# _ # _
# |\/| / \ _ # |\/| / \ _
@ -131,24 +139,20 @@ print "END_BASIS_SET"
# #
# Function # Function
# #
d_gms_order ={ 0:["s"], d_gms_order = ["s",
1:[ "x", "y", "z" ], "x", "y", "z" ,
2:[ "xx", "yy", "zz", "xy", "xz", "yz" ], "xx", "yy", "zz", "xy", "xz", "yz" ,
3:[ "xxx", "yyy", "zzz", "xxy", "xxz", "yyx", "yyz", "zzx", "zzy", "xyz"], "xxx", "yyy", "zzz", "xxy", "xxz", "yyx", "yyz", "zzx", "zzy", "xyz",
4:[ "xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz", "zzxy"] } "xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz", "zzxy"]
def compare_gamess_style(item1, item2): def key_gamess_style(item1):
n1,n2 = map(len,(item1,item2))
assert (n1 == n2)
try: try:
l = d_gms_order[n1] l = d_gms_order.index(item1)
except KeyError: except ValueError:
return 0 return 0
# raise (KeyError, "We dont handle L than 4")
else: else:
a = l.index(item1) result = d_gms_order.index(item1)
b = l.index(item2) return result
return cmp( a, b )
def expend_sym_str(str_): def expend_sym_str(str_):
#Expend x2 -> xx #Expend x2 -> xx
@ -184,9 +188,10 @@ def get_nb_permutation(str_):
if (str_) == 's': return 1 if (str_) == 's': return 1
else: return n_orbital(len(str_)) else: return n_orbital(len(str_))
import functools
def order_l_l_sym(l_l_sym): def order_l_l_sym(l_l_sym):
n = 1 n = 1
iter_ = range(len(l_l_sym)) iter_ = list(range(len(l_l_sym)))
for i in iter_: for i in iter_:
if n != 1: if n != 1:
n += -1 n += -1
@ -195,9 +200,11 @@ def order_l_l_sym(l_l_sym):
l = l_l_sym[i] l = l_l_sym[i]
n = get_nb_permutation(l[2]) n = get_nb_permutation(l[2])
def local_key(x):
return key_gamess_style(x[2])
l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n], l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n],
key=lambda x: x[2], key=local_key)
cmp=compare_gamess_style)
return l_l_sym return l_l_sym
@ -237,21 +244,21 @@ def order_phase(mo_coef):
def chunked(l, chunks_size): def chunked(l, chunks_size):
l_block = [] l_block = []
for i in l: for i in l:
chunks = [i[x:x + chunks_size] for x in xrange(0, len(i), chunks_size)] chunks = [i[x:x + chunks_size] for x in range(0, len(i), chunks_size)]
l_block.append(chunks) l_block.append(chunks)
return l_block return l_block
def print_mo_coef(mo_coef_block, l_l_sym): def print_mo_coef(mo_coef_block, l_l_sym):
print "" print("")
print "BEGIN_MO" print("BEGIN_MO")
print "" print("")
len_block_curent = 0 len_block_curent = 0
nb_block = len(mo_coef_block[0]) nb_block = len(mo_coef_block[0])
for i_block in range(0, nb_block): for i_block in range(0, nb_block):
a = [i[i_block] for i in mo_coef_block] a = [i[i_block] for i in mo_coef_block]
r_ = range(len_block_curent, len_block_curent + len(a[0])) r_ = list(range(len_block_curent, len_block_curent + len(a[0])))
print " ".join([str(i + 1) for i in r_]) print(" ".join([str(i + 1) for i in r_]))
len_block_curent += len(a[0]) len_block_curent += len(a[0])
@ -260,17 +267,17 @@ def print_mo_coef(mo_coef_block, l_l_sym):
i_a = int(l[1]) - 1 i_a = int(l[1]) - 1
sym = l[2] sym = l[2]
print l_label[i_a], sym, " ".join('%20.15e'%i print(l_label[i_a], sym, " ".join('%20.15e'%i
for i in a[i]) for i in a[i]))
if i_block != nb_block - 1: if i_block != nb_block - 1:
print "" print("")
else: else:
print "END_MO" print("END_MO")
mo_coef = ezfio.get_mo_basis_mo_coef() mo_coef = ezfio.get_mo_basis_mo_coef()
mo_coef_transp = zip(*mo_coef) mo_coef_transp = list(zip(*mo_coef))
mo_coef_block = chunked(mo_coef_transp, 4) mo_coef_block = chunked(mo_coef_transp, 4)
print_mo_coef(mo_coef_block, l_l_sym_ordered) print_mo_coef(mo_coef_block, l_l_sym_ordered)
@ -279,8 +286,8 @@ print_mo_coef(mo_coef_block, l_l_sym_ordered)
# | _> (/_ |_| (_| (_) # | _> (/_ |_| (_| (_)
# #
if do_pseudo: if do_pseudo:
print "" print("")
print "BEGIN_PSEUDO" print("BEGIN_PSEUDO")
klocmax = ezfio.get_pseudo_pseudo_klocmax() klocmax = ezfio.get_pseudo_pseudo_klocmax()
kmax = ezfio.get_pseudo_pseudo_kmax() kmax = ezfio.get_pseudo_pseudo_kmax()
lmax = ezfio.get_pseudo_pseudo_lmax() lmax = ezfio.get_pseudo_pseudo_lmax()
@ -318,17 +325,17 @@ if do_pseudo:
l_str.append(l_dump) l_str.append(l_dump)
str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE" str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE"
print str_.format(a, i + 1, int(d_z[a])-int(l_charge[i]), int(len(l_str) - 1)) print(str_.format(a, i + 1, int(d_z[a])-int(l_charge[i]), int(len(l_str) - 1)))
for i, l in enumerate(l_str): for i, l in enumerate(l_str):
str_ = "FOR L= {0} COEFF N ZETA" str_ = "FOR L= {0} COEFF N ZETA"
print str_.format(int(len(l_str) - i - 1)) print(str_.format(int(len(l_str) - i - 1)))
for ii, ll in enumerate(l): for ii, ll in enumerate(l):
print " ", ii + 1, ll print(" ", ii + 1, ll)
str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS." str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS."
print str_.format(sum([int(d_z[a])-int(l_charge[i]) for i,a in enumerate(l_label)])) print(str_.format(sum([int(d_z[a])-int(l_charge[i]) for i,a in enumerate(l_label)])))
print "END_PSEUDO" print("END_PSEUDO")
# _ # _
# | \ _ _|_ # | \ _ _|_
@ -343,12 +350,12 @@ bit_kind = ezfio.get_determinants_bit_kind()
nexcitedstate = ezfio.get_determinants_n_states() nexcitedstate = ezfio.get_determinants_n_states()
print "" print("")
print "BEGIN_DET" print("BEGIN_DET")
print "" print("")
print "mo_num", mo_num print("mo_num", mo_num)
print "det_num", n_det print("det_num", n_det)
print "" print("")
if "QP_STATE" in os.environ: if "QP_STATE" in os.environ:
state = int(os.environ["QP_STATE"])-1 state = int(os.environ["QP_STATE"])-1
@ -372,13 +379,13 @@ MultiDetAlpha = []
MultiDetBeta = [] MultiDetBeta = []
for coef, (det_a, det_b) in zip(psi_coef_small, psi_det): for coef, (det_a, det_b) in zip(psi_coef_small, psi_det):
print coef print(coef)
MyDetA=decode(det_a) MyDetA=decode(det_a)
MyDetB=decode(det_b) MyDetB=decode(det_b)
print MyDetA print(MyDetA)
print MyDetB print(MyDetB)
print '' print('')
MultiDetAlpha.append( det_a ) MultiDetAlpha.append( det_a )
MultiDetBeta.append( det_b ) MultiDetBeta.append( det_b )
print "END_DET" print("END_DET")

View File

@ -23,73 +23,73 @@ program qmcpack
call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5') call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5')
call system('$QP_ROOT/src/champ/qp_convert.py '//trim(ezfio_filename)) call system('$QP_ROOT/src/champ/qp_convert.py '//trim(ezfio_filename))
integer :: iunit ! integer :: iunit
integer, external :: getUnitAndOpen ! integer, external :: getUnitAndOpen
iunit = getUnitAndOpen(trim(ezfio_filename)//'.H','w') ! iunit = getUnitAndOpen(trim(ezfio_filename)//'.H','w')
!
double precision, external :: diag_h_mat_elem ! double precision, external :: diag_h_mat_elem
write(iunit,*) N_states ! write(iunit,*) N_states
do istate=1,N_states ! do istate=1,N_states
write(iunit,*) istate, psi_energy_with_nucl_rep(istate) ! write(iunit,*) istate, psi_energy_with_nucl_rep(istate)
enddo ! enddo
write(iunit,*) N_det ! write(iunit,*) N_det
do k=1,N_det ! do k=1,N_det
write(iunit,'(I10,X,F22.15)') k, diag_h_mat_elem(psi_det(1,1,k),N_int) + nuclear_repulsion ! write(iunit,'(I10,X,F22.15)') k, diag_h_mat_elem(psi_det(1,1,k),N_int) + nuclear_repulsion
enddo ! enddo
!
!
double precision :: F(N_states) ! double precision :: F(N_states)
integer(bit_kind), allocatable :: det(:,:,:) ! integer(bit_kind), allocatable :: det(:,:,:)
double precision , allocatable :: coef(:,:) ! double precision , allocatable :: coef(:,:)
integer :: ispin ! integer :: ispin
double precision :: norm(N_states), hij ! double precision :: norm(N_states), hij
allocate(det(N_int,2,N_det), coef(N_det,N_states)) ! allocate(det(N_int,2,N_det), coef(N_det,N_states))
do j=1,mo_num ! do j=1,mo_num
do i=1,j-1 ! do i=1,j-1
do ispin=1,2 ! do ispin=1,2
call build_singly_excited_wavefunction(j,i,1,det,coef) ! call build_singly_excited_wavefunction(j,i,1,det,coef)
F = 0.d0 ! F = 0.d0
do istate=1,N_states ! do istate=1,N_states
norm(istate) = 0.d0 ! norm(istate) = 0.d0
do k=1,N_det ! do k=1,N_det
norm(istate) = norm(istate) + coef(k,istate) * coef(k,istate) ! norm(istate) = norm(istate) + coef(k,istate) * coef(k,istate)
enddo ! enddo
if (norm(istate) > 0.d0) then ! if (norm(istate) > 0.d0) then
norm(istate) = (1.d0/dsqrt(norm(istate))) ! norm(istate) = (1.d0/dsqrt(norm(istate)))
endif ! endif
enddo ! enddo
if (sum(norm(:)) > 0.d0) then ! if (sum(norm(:)) > 0.d0) then
!
do istate = 1,N_states ! do istate = 1,N_states
coef(:,istate) = coef(:,istate) * norm(istate) ! coef(:,istate) = coef(:,istate) * norm(istate)
enddo ! enddo
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,istate,hij) REDUCTION(+:F) ! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,istate,hij) REDUCTION(+:F)
do k=1,N_det ! do k=1,N_det
if (sum(coef(k,:)*coef(k,:))==0.d0) cycle ! if (sum(coef(k,:)*coef(k,:))==0.d0) cycle
call i_H_j(det(1,1,k), det(1,1,k), N_int, hij) ! call i_H_j(det(1,1,k), det(1,1,k), N_int, hij)
do istate=1,N_states ! do istate=1,N_states
F(istate) = F(istate) + hij*coef(k,istate)*coef(k,istate) ! F(istate) = F(istate) + hij*coef(k,istate)*coef(k,istate)
enddo ! enddo
do l=1,k-1 ! do l=1,k-1
if (sum(coef(l,:)*coef(l,:))==0.d0) cycle ! if (sum(coef(l,:)*coef(l,:))==0.d0) cycle
call i_H_j(det(1,1,k), det(1,1,l), N_int, hij) ! call i_H_j(det(1,1,k), det(1,1,l), N_int, hij)
do istate=1,N_states ! do istate=1,N_states
F(istate) = F(istate) + 2.d0*hij*coef(k,istate)*coef(l,istate) ! F(istate) = F(istate) + 2.d0*hij*coef(k,istate)*coef(l,istate)
enddo ! enddo
enddo ! enddo
enddo ! enddo
!$OMP END PARALLEL DO ! !$OMP END PARALLEL DO
F(:) = F(:) - psi_energy(:) ! F(:) = F(:) - psi_energy(:)
endif ! endif
do istate=1,N_states ! do istate=1,N_states
write(iunit,'(I4,X,I4,X,I1,X,I3,X,F22.15)') i, j, ispin, istate, F(istate) ! write(iunit,'(I4,X,I4,X,I1,X,I3,X,F22.15)') i, j, ispin, istate, F(istate)
enddo ! enddo
enddo ! enddo
enddo ! enddo
enddo ! enddo
!
deallocate(det,coef) ! deallocate(det,coef)
!
close(iunit) ! close(iunit)
end end