PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-14 10:03:47 +01:00
Code Releases Activity

Merge branch 'QuantumPackage-dev-stable' into dev-stable-tc-scf

Browse Source
This commit is contained in:
Abdallah Ammar 2023-05-06 20:42:52 +02:00
commit b74ac64148
54 changed files with 1017 additions and 620 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort_2021_xHost.cfg
View File

@ -6,7 +6,7 @@
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
FC : ifort -fpic -diag-disable 5462
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=64 -DINTEL

2
configure vendored
View File

@ -232,7 +232,7 @@ EOF
execute << EOF
cd "\${QP_ROOT}"/external
tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.2.tar.gz
tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.?.tar.gz
cd f77-zmq-*
./configure --prefix=\$QP_ROOT
export ZMQ_H="\$QP_ROOT"/include/zmq.h

2
scripts/compilation/qp_create_ninja
View File

@ -25,7 +25,7 @@ except ImportError:
"quantum_package.rc"))
print("\n".join(["", "Error:", "source %s" % f, ""]))
sys.exit(1)
raise
# Compress path
def comp_path(path):

186
scripts/qp_exc_energy.py Executable file
View File

@ -0,0 +1,186 @@
#!/usr/bin/env python
# Computes the error on the excitation energy of a CIPSI run.
def student(p,df):
import scipy
from scipy.stats import t
return t.ppf(p, df)
def chi2cdf(x,k):
import scipy
import scipy.stats
return scipy.stats.chi2.cdf(x,k)
def jarque_bera(data):
n = max(len(data), 2)
norm = 1./ sum( [ w for (_,w) in data ] )
mu = sum( [ w* x for (x,w) in data ] ) * norm
sigma2 = sum( [ w*(x-mu)**2 for (x,w) in data ] ) * norm
if sigma2 > 0.:
S = ( sum( [ w*(x-mu)**3 for (x,w) in data ] ) * norm ) / sigma2**(3./2.)
K = ( sum( [ w*(x-mu)**4 for (x,w) in data ] ) * norm ) / sigma2**2
else:
S = 0.
K = 0.
# Value of the Jarque-Bera test
JB = n/6. * (S**2 + 1./4. * (K-3.)**2)
# Probability that the data comes from a Gaussian distribution
p = 1. - chi2cdf(JB,2)
return JB, mu, sqrt(sigma2/(n-1)), p
to_eV = 27.2107362681
import sys, os
import scipy
import scipy.stats
from math import sqrt, gamma, exp
import json
def read_data(filename,state):
""" Read energies and PT2 from input file """
with open(filename,'r') as f:
lines = json.load(f)['fci']
print(f"State: {state}")
gs = []
es = []
for l in lines:
try:
pt2_0 = l['states'][0]['pt2']
e_0 = l['states'][0]['energy']
pt2_1 = l['states'][state]['pt2']
e_1 = l['states'][state]['energy']
gs.append( (e_0, pt2_0) )
es.append( (e_1, pt2_1) )
except: pass
def f(p_1, p0, p1):
e, pt2 = p0
y0, x0 = p_1
y1, x1 = p1
try:
alpha = (y1-y0)/(x0-x1)
except ZeroDivisionError:
alpha = 1.
return [e, pt2, alpha]
for l in (gs, es):
p_1, p0, p1 = l[0], l[0], l[1]
l[0] = f(p_1, p0, p1)
for i in range(1,len(l)-1):
p_1 = (l[i-1][0], l[i-1][1])
p0 = l[i]
p1 = l[i+1]
l[i] = f(p_1, p0, p1)
i = len(l)-1
p_1 = (l[i-1][0], l[i-1][1])
p0 = l[i]
p1 = l[-1]
l[i] = f(p_1, p0, p1)
return [ x+y for x,y in zip(gs,es) ]
def compute(data):
d = []
for e0, p0, a0, e1, p1, a1 in data:
x = (e1+p1)-(e0+p0)
w = 1./sqrt(p0**2 + p1**2)
bias = (a1-1.)*p1 - (a0-1.)*p0
d.append( (x,w,bias) )
x = []
target = (scipy.stats.norm.cdf(1.)-0.5)*2 # = 0.6827
print("| %2s | %8s | %8s | %8s | %8s | %8s |"%( "N", "DE", "+/-", "bias", "P(G)", "J"))
print("|----+----------+----------+----------+----------+----------|")
xmax = (0.,0.,0.,0.,0.,0,0.)
for i in range(len(data)-1):
jb, mu, sigma, p = jarque_bera( [ (x,w) for (x,w,bias) in d[i:] ] )
bias = sum ( [ w * e for (_,w,e) in d[i:] ] ) / sum ( [ w for (_,w,_) in d[i:] ] )
mu = (mu+0.5*bias) * to_eV
sigma = sigma * to_eV
bias = bias * to_eV
n = len(data[i:])
beta = student(0.5*(1.+target/p) ,n)
err = sigma * beta + 0.5*abs(bias)
print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |"%( n, mu, err, bias, p, jb))
if n < 3 :
continue
y = (err, p, mu, err, jb,n,bias)
if p > xmax[1]: xmax = y
if p < 0.8:
continue
x.append(y)
x = sorted(x)
print("|----+----------+----------+----------+----------+----------|")
if x != []:
xmax = x[0]
_, p, mu, err, jb, n, bias = xmax
beta = student(0.5*(1.+target/p),n)
print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |\n"%(n, mu, err, bias, p, jb))
return mu, err, bias, p
filename = sys.argv[1]
print(filename)
if len(sys.argv) > 2:
state = int(sys.argv[2])
else:
state = 1
data = read_data(filename,state)
mu, err, bias, _ = compute(data)
print(" %s: %8.3f +/- %5.3f eV\n"%(filename, mu, err))
import numpy as np
A = np.array( [ [ data[-1][1], 1. ],
[ data[-2][1], 1. ] ] )
B = np.array( [ [ data[-1][0] ],
[ data[-2][0] ] ] )
E0 = np.linalg.solve(A,B)[1]
A = np.array( [ [ data[-1][4], 1. ],
[ data[-2][4], 1. ] ] )
B = np.array( [ [ data[-1][3] ],
[ data[-2][3] ] ] )
E1 = np.linalg.solve(A,B)[1]
average_2 = (E1-E0)*to_eV
A = np.array( [ [ data[-1][1], 1. ],
[ data[-2][1], 1. ],
[ data[-3][1], 1. ] ] )
B = np.array( [ [ data[-1][0] ],
[ data[-2][0] ],
[ data[-3][0] ] ] )
E0 = np.linalg.lstsq(A,B,rcond=None)[0][1]
A = np.array( [ [ data[-1][4], 1. ],
[ data[-2][4], 1. ],
[ data[-3][4], 1. ] ] )
B = np.array( [ [ data[-1][3] ],
[ data[-2][3] ],
[ data[-3][3] ] ] )
E1 = np.linalg.lstsq(A,B,rcond=None)[0][1]
average_3 = (E1-E0)*to_eV
exc = ((data[-1][3] + data[-1][4]) - (data[-1][0] + data[-1][1])) * to_eV
error_2 = abs(average_2 - average_3)
error_3 = abs(average_3 - exc)
print(" 2-3 points: %.3f +/- %.3f "% (average_3, error_2))
print(" largest wf: %.3f +/- %.3f "%(average_3, error_3))

5
src/ao_two_e_ints/EZFIO.cfg
View File

@ -18,3 +18,8 @@ interface: ezfio,provider,ocaml
default: False
ezfio_name: direct
[do_ao_cholesky]
type: logical
doc: Perform Cholesky decomposition of AO integrals
interface: ezfio,provider,ocaml
default: True

100
src/ao_two_e_ints/cholesky.irp.f Normal file
View File

@ -0,0 +1,100 @@
BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
implicit none
BEGIN_DOC
! Number of Cholesky vectors in AO basis
END_DOC
integer :: i,j,k,l
double precision :: xnorm0, x, integral
double precision, external :: ao_two_e_integral
cholesky_ao_num_guess = 0
xnorm0 = 0.d0
x = 0.d0
do j=1,ao_num
do i=1,ao_num
integral = ao_two_e_integral(i,i,j,j)
if (integral > ao_integrals_threshold) then
cholesky_ao_num_guess += 1
else
x += integral
endif
enddo
enddo
print *, 'Cholesky decomposition of AO integrals'
print *, '--------------------------------------'
print *, ''
print *, 'Estimated Error: ', x
print *, 'Guess size: ', cholesky_ao_num_guess, '(', 100.d0*dble(cholesky_ao_num_guess)/dble(ao_num*ao_num), ' %)'
END_PROVIDER
BEGIN_PROVIDER [ integer, cholesky_ao_num ]
&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, cholesky_ao_num_guess) ]
use mmap_module
implicit none
BEGIN_DOC
! Cholesky vectors in AO basis: (ik|a):
! <ij|kl> = (ik|jl) = sum_a (ik|a).(a|jl)
END_DOC
type(c_ptr) :: ptr
integer :: fd, i,j,k,l, rank
double precision, pointer :: ao_integrals(:,:,:,:)
double precision, external :: ao_two_e_integral
! Store AO integrals in a memory mapped file
call mmap(trim(ezfio_work_dir)//'ao_integrals', &
(/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), &
8, fd, .False., ptr)
call c_f_pointer(ptr, ao_integrals, (/ao_num, ao_num, ao_num, ao_num/))
double precision :: integral
logical, external :: ao_two_e_integral_zero
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l, integral) SCHEDULE(dynamic)
do l=1,ao_num
do j=1,l
do k=1,ao_num
do i=1,k
if (ao_two_e_integral_zero(i,j,k,l)) cycle
integral = ao_two_e_integral(i,k,j,l)
ao_integrals(i,k,j,l) = integral
ao_integrals(k,i,j,l) = integral
ao_integrals(i,k,l,j) = integral
ao_integrals(k,i,l,j) = integral
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
! Call Lapack
cholesky_ao_num = cholesky_ao_num_guess
call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_integrals_threshold, ao_num*ao_num, cholesky_ao)
print *, 'Rank: ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
! Remove mmap
double precision, external :: getUnitAndOpen
call munmap( &
(/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), &
8, fd, ptr)
open(unit=99,file=trim(ezfio_work_dir)//'ao_integrals')
close(99, status='delete')
END_PROVIDER
BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num, ao_num) ]
implicit none
BEGIN_DOC
! Transposed of the Cholesky vectors in AO basis set
END_DOC
integer :: i,j,k
do j=1,ao_num
do i=1,ao_num
do k=1,ao_num
cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k)
enddo
enddo
enddo
END_PROVIDER

2
src/ao_two_e_ints/map_integrals.irp.f
View File

@ -486,7 +486,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
PROVIDE ao_two_e_integrals_in_map ao_integrals_map
if (ao_one_e_integral_zero(j,l)) then
out_val = 0.d0
out_val(1:sze) = 0.d0
return
endif

1
src/cipsi/NEED
View File

@ -1,3 +1,4 @@
json
perturbation
zmq
mpi

24
src/cipsi/cipsi.irp.f
View File

@ -16,7 +16,6 @@ subroutine run_cipsi
double precision, external :: memory_of_double
PROVIDE H_apply_buffer_allocated
N_iter = 1
threshold_generators = 1.d0
SOFT_TOUCH threshold_generators
@ -76,7 +75,6 @@ subroutine run_cipsi
)
write(*,'(A)') '--------------------------------------------------------------------------------'
to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
to_select = max(N_states_diag, to_select)
if (do_pt2) then
@ -106,10 +104,10 @@ subroutine run_cipsi
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
call print_extrapolated_energy()
call print_mol_properties()
N_iter += 1
call write_cipsi_json(pt2_data,pt2_data_err)
if (qp_stop()) exit
@ -130,13 +128,13 @@ subroutine run_cipsi
if (qp_stop()) exit
enddo
if (.not.qp_stop()) then
if (N_det < N_det_max) then
call diagonalize_CI
call save_wavefunction
call save_energy(psi_energy_with_nucl_rep, zeros)
endif
! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
if ((.not.qp_stop()).and. &
(N_det > N_det_max) .and. &
(maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. &
(maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
(correlation_energy_ratio <= correlation_energy_ratio_max) &
) then
if (do_pt2) then
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
@ -155,11 +153,13 @@ subroutine run_cipsi
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
call print_summary(psi_energy_with_nucl_rep(1:N_states), &
pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
call print_extrapolated_energy()
call print_mol_properties()
call write_cipsi_json(pt2_data,pt2_data_err)
endif
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
end

2
src/cipsi/environment.irp.f
View File

@ -7,7 +7,9 @@ BEGIN_PROVIDER [ integer, nthreads_pt2 ]
character*(32) :: env
call getenv('QP_NTHREADS_PT2',env)
if (trim(env) /= '') then
call lock_io()
read(env,*) nthreads_pt2
call unlock_io()
call write_int(6,nthreads_pt2,'Target number of threads for PT2')
endif
END_PROVIDER

22
src/cipsi/stochastic_cipsi.irp.f
View File

@ -15,7 +15,6 @@ subroutine run_stochastic_cipsi
double precision, external :: memory_of_double
PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map
N_iter = 1
threshold_generators = 1.d0
SOFT_TOUCH threshold_generators
@ -96,10 +95,10 @@ subroutine run_stochastic_cipsi
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
call print_extrapolated_energy()
call print_mol_properties()
N_iter += 1
call write_cipsi_json(pt2_data,pt2_data_err)
if (qp_stop()) exit
@ -119,13 +118,13 @@ subroutine run_stochastic_cipsi
if (qp_stop()) exit
enddo
if (.not.qp_stop()) then
if (N_det < N_det_max) then
call diagonalize_CI
call save_wavefunction
call save_energy(psi_energy_with_nucl_rep, zeros)
endif
! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
if ((.not.qp_stop()).and. &
(N_det > N_det_max) .and. &
(maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. &
(maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
(correlation_energy_ratio <= correlation_energy_ratio_max) &
) then
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states)
@ -135,9 +134,10 @@ subroutine run_stochastic_cipsi
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
call print_summary(psi_energy_with_nucl_rep, &
pt2_data , pt2_data_err, N_det, N_configuration, N_states, psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
call print_extrapolated_energy()
call print_mol_properties()
call write_cipsi_json(pt2_data,pt2_data_err)
endif
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)

53
src/cipsi/write_cipsi_json.irp.f Normal file
View File

@ -0,0 +1,53 @@
subroutine write_cipsi_json(pt2_data, pt2_data_err)
use selection_types
implicit none
BEGIN_DOC
! Writes JSON data for CIPSI runs
END_DOC
type(pt2_type), intent(in) :: pt2_data, pt2_data_err
integer :: i,j,k
call lock_io
character*(64), allocatable :: fmtk(:)
integer :: N_states_p, N_iter_p
N_states_p = min(N_states,N_det)
N_iter_p = min(N_iter,8)
allocate(fmtk(0:N_iter_p))
fmtk(:) = '('' '',E22.15,'','')'
fmtk(N_iter_p) = '('' '',E22.15)'
write(json_unit, json_dict_uopen_fmt)
write(json_unit, json_int_fmt) 'n_det', N_det
if (s2_eig) then
write(json_unit, json_int_fmt) 'n_cfg', N_configuration
if (only_expected_s2) then
write(json_unit, json_int_fmt) 'n_csf', N_csf
endif
endif
write(json_unit, json_array_open_fmt) 'states'
do k=1,N_states_p
write(json_unit, json_dict_uopen_fmt)
write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k)
write(json_unit, json_real_fmt) 's2', psi_s2(k)
write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k)
write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k)
write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k)
write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k)
write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k)
write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k)
write(json_unit, json_array_open_fmt) 'ex_energy'
do i=2,N_iter_p
write(json_unit, fmtk(i)) extrapolated_energy(i,k)
enddo
write(json_unit, json_array_close_fmtx)
if (k < N_states_p) then
write(json_unit, json_dict_close_fmt)
else
write(json_unit, json_dict_close_fmtx)
endif
enddo
write(json_unit, json_array_close_fmtx)
write(json_unit, json_dict_close_fmt)
deallocate(fmtk)
call unlock_io
end

3
src/cipsi_tc_bi_ortho/NEED
View File

@ -1,6 +1,7 @@
json
mpi
perturbation
zmq
iterations_tc
iterations
csf
tc_bi_ortho

10
src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
View File

@ -94,7 +94,15 @@ subroutine run_stochastic_cipsi
call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
! stop
N_iter += 1
call print_summary(psi_energy_with_nucl_rep, &
pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
call print_extrapolated_energy()
! call print_mol_properties()
call write_cipsi_json(pt2_data,pt2_data_err)
if (qp_stop()) exit

53
src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f Normal file
View File

@ -0,0 +1,53 @@
subroutine write_cipsi_json(pt2_data, pt2_data_err)
use selection_types
implicit none
BEGIN_DOC
! Writes JSON data for CIPSI runs
END_DOC
type(pt2_type), intent(in) :: pt2_data, pt2_data_err
integer :: i,j,k
call lock_io
character*(64), allocatable :: fmtk(:)
integer :: N_states_p, N_iter_p
N_states_p = min(N_states,N_det)
N_iter_p = min(N_iter,8)
allocate(fmtk(0:N_iter_p))
fmtk(:) = '('' '',E22.15,'','')'
fmtk(N_iter_p) = '('' '',E22.15)'
write(json_unit, json_dict_uopen_fmt)
write(json_unit, json_int_fmt) 'n_det', N_det
if (s2_eig) then
write(json_unit, json_int_fmt) 'n_cfg', N_configuration
if (only_expected_s2) then
write(json_unit, json_int_fmt) 'n_csf', N_csf
endif
endif
write(json_unit, json_array_open_fmt) 'states'
do k=1,N_states_p
write(json_unit, json_dict_uopen_fmt)
write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k)
write(json_unit, json_real_fmt) 's2', psi_s2(k)
write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k)
write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k)
write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k)
write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k)
write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k)
write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k)
write(json_unit, json_array_open_fmt) 'ex_energy'
do i=2,N_iter_p
write(json_unit, fmtk(i)) extrapolated_energy(i,k)
enddo
write(json_unit, json_array_close_fmtx)
if (k < N_states_p) then
write(json_unit, json_dict_close_fmt)
else
write(json_unit, json_dict_close_fmtx)
endif
enddo
write(json_unit, json_array_close_fmtx)
write(json_unit, json_dict_close_fmt)
deallocate(fmtk)
call unlock_io
end

4
src/davidson/davidson_parallel.irp.f
View File

@ -150,7 +150,9 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze,
exit
endif
if(task_id == 0) exit
call lock_io()
read (msg,*) imin, imax, ishift, istep
call unlock_io()
integer :: k
do k=imin,imax
v_t(:,k) = 0.d0
@ -546,6 +548,8 @@ end
integer function zmq_put_N_states_diag(zmq_to_qp_run_socket,worker_id)
use f77_zmq
implicit none

3
src/davidson/diagonalization_hs2_dressed.irp.f
View File

@ -461,9 +461,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
integer :: lwork, info
double precision, allocatable :: work(:)
y = h
!y = h_p
! y = h_p ! Doesn't work for non-singlets
lwork = -1
allocate(work(1))
call dsygv(1,'V','U',shift2,y,size(y,1), &

17
src/ezfio_files/lock.irp.f
View File

@ -9,4 +9,21 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), file_lock ]
call omp_init_lock(file_lock)
END_PROVIDER
! These functions need to be called because internal read and write are not thread safe.
subroutine lock_io()
implicit none
BEGIN_DOC
! Needs to be called because before doing I/O because internal read and write
! are not thread safe.
END_DOC
call omp_set_lock(file_lock)
end subroutine lock_io()
subroutine unlock_io()
implicit none
BEGIN_DOC
! Needs to be called because afterdoing I/O because internal read and write
! are not thread safe.
END_DOC
call omp_unset_lock(file_lock)
end subroutine lock_io()

7
src/fci/fci.irp.f
View File

@ -39,12 +39,19 @@ program fci
if (.not.is_zmq_slave) then
PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
write(json_unit,json_array_open_fmt) 'fci'
if (do_pt2) then
call run_stochastic_cipsi
else
call run_cipsi
endif
write(json_unit,json_dict_uopen_fmt)
write(json_unit,json_dict_close_fmtx)
write(json_unit,json_array_close_fmtx)
call json_close
else
PROVIDE mo_two_e_integrals_in_map pt2_min_parallel_tasks

1
src/fci_tc_bi/NEED
View File

@ -1,3 +1,4 @@
json
tc_bi_ortho
davidson_undressed
cipsi_tc_bi_ortho

6
src/fci_tc_bi/fci_tc_bi_ortho.irp.f
View File

@ -62,6 +62,7 @@ subroutine run_cipsi_tc
endif
endif
! ---
write(json_unit,json_array_open_fmt) 'fci_tc'
if (do_pt2) then
call run_stochastic_cipsi
@ -69,6 +70,11 @@ subroutine run_cipsi_tc
call run_cipsi
endif
write(json_unit,json_dict_uopen_fmt)
write(json_unit,json_dict_close_fmtx)
write(json_unit,json_array_close_fmtx)
call json_close
else
PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks
if(elec_alpha_num+elec_beta_num.ge.3)then

4
src/fci_tc_bi/save_energy.irp.f
View File

@ -4,6 +4,6 @@ subroutine save_energy(E,pt2)
! Saves the energy in |EZFIO|.
END_DOC
double precision, intent(in) :: E(N_states), pt2(N_states)
call ezfio_set_fci_tc_energy(E(1:N_states))
call ezfio_set_fci_tc_energy_pt2(E(1:N_states)+pt2(1:N_states))
call ezfio_set_fci_tc_bi_energy(E(1:N_states))
call ezfio_set_fci_tc_bi_energy_pt2(E(1:N_states)+pt2(1:N_states))
end

1
src/hartree_fock/NEED
View File

@ -1,3 +1,4 @@
ao_one_e_ints
ao_two_e_ints
scf_utils
json

311
src/hartree_fock/fock_matrix_hf.irp.f
View File

@ -15,115 +15,59 @@
double precision, allocatable :: ao_two_e_integral_alpha_tmp(:,:)
double precision, allocatable :: ao_two_e_integral_beta_tmp(:,:)
ao_two_e_integral_alpha = 0.d0
ao_two_e_integral_beta = 0.d0
if (do_direct_integrals) then
if (do_ao_cholesky) then ! Use Cholesky-decomposed integrals
ao_two_e_integral_alpha(:,:) = ao_two_e_integral_alpha_chol(:,:)
ao_two_e_integral_beta (:,:) = ao_two_e_integral_beta_chol (:,:)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0, &
!$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2, &
!$OMP local_threshold)&
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
!$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, &
!$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta)
else ! Use integrals in AO basis set
allocate(keys(1), values(1))
allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), &
ao_two_e_integral_beta_tmp(ao_num,ao_num))
ao_two_e_integral_alpha_tmp = 0.d0
ao_two_e_integral_beta_tmp = 0.d0
ao_two_e_integral_alpha = 0.d0
ao_two_e_integral_beta = 0.d0
if (do_direct_integrals) then
q = ao_num*ao_num*ao_num*ao_num
!$OMP DO SCHEDULE(static,64)
do p=1_8,q
call two_e_integrals_index_reverse(kk,ii,ll,jj,p)
if ( (kk(1)>ao_num).or. &
(ii(1)>ao_num).or. &
(jj(1)>ao_num).or. &
(ll(1)>ao_num) ) then
cycle
endif
k = kk(1)
i = ii(1)
l = ll(1)
j = jj(1)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0,&
!$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2,&
!$OMP local_threshold) &
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
!$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz,&
!$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta)
logical, external :: ao_two_e_integral_zero
if (ao_two_e_integral_zero(i,k,j,l)) then
cycle
endif
local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j)
if (local_threshold < ao_integrals_threshold) then
cycle
endif
i0 = i
j0 = j
k0 = k
l0 = l
values(1) = 0.d0
local_threshold = ao_integrals_threshold/local_threshold
do k2=1,8
if (kk(k2)==0) then
cycle
endif
i = ii(k2)
j = jj(k2)
k = kk(k2)
l = ll(k2)
c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)
c1 = SCF_density_matrix_ao_alpha(k,i)
c2 = SCF_density_matrix_ao_beta(k,i)
if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then
cycle
endif
if (values(1) == 0.d0) then
values(1) = ao_two_e_integral(k0,l0,i0,j0)
endif
integral = c0 * values(1)
ao_two_e_integral_alpha_tmp(i,j) += integral
ao_two_e_integral_beta_tmp (i,j) += integral
integral = values(1)
ao_two_e_integral_alpha_tmp(l,j) -= c1 * integral
ao_two_e_integral_beta_tmp (l,j) -= c2 * integral
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp
ao_two_e_integral_beta += ao_two_e_integral_beta_tmp
!$OMP END CRITICAL
deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)
!$OMP END PARALLEL
else
PROVIDE ao_two_e_integrals_in_map
allocate(keys(1), values(1))
allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), &
ao_two_e_integral_beta_tmp(ao_num,ao_num))
ao_two_e_integral_alpha_tmp = 0.d0
ao_two_e_integral_beta_tmp = 0.d0
integer(omp_lock_kind) :: lck(ao_num)
integer(map_size_kind) :: i8
integer :: ii(8), jj(8), kk(8), ll(8), k2
integer(cache_map_size_kind) :: n_elements_max, n_elements
integer(key_kind), allocatable :: keys(:)
double precision, allocatable :: values(:)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, &
!$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)&
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
!$OMP ao_integrals_map, ao_two_e_integral_alpha, ao_two_e_integral_beta)
call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max)
allocate(keys(n_elements_max), values(n_elements_max))
allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), &
ao_two_e_integral_beta_tmp(ao_num,ao_num))
ao_two_e_integral_alpha_tmp = 0.d0
ao_two_e_integral_beta_tmp = 0.d0
!$OMP DO SCHEDULE(static,1)
do i8=0_8,ao_integrals_map%map_size
n_elements = n_elements_max
call get_cache_map(ao_integrals_map,i8,keys,values,n_elements)
do k1=1,n_elements
call two_e_integrals_index_reverse(kk,ii,ll,jj,keys(k1))
q = ao_num*ao_num*ao_num*ao_num
!$OMP DO SCHEDULE(static,64)
do p=1_8,q
call two_e_integrals_index_reverse(kk,ii,ll,jj,p)
if ( (kk(1)>ao_num).or. &
(ii(1)>ao_num).or. &
(jj(1)>ao_num).or. &
(ll(1)>ao_num) ) then
cycle
endif
k = kk(1)
i = ii(1)
l = ll(1)
j = jj(1)
logical, external :: ao_two_e_integral_zero
if (ao_two_e_integral_zero(i,k,j,l)) then
cycle
endif
local_threshold = ao_two_e_integral_schwartz(k,l)*ao_two_e_integral_schwartz(i,j)
if (local_threshold < ao_integrals_threshold) then
cycle
endif
i0 = i
j0 = j
k0 = k
l0 = l
values(1) = 0.d0
local_threshold = ao_integrals_threshold/local_threshold
do k2=1,8
if (kk(k2)==0) then
cycle
@ -132,25 +76,162 @@
j = jj(k2)
k = kk(k2)
l = ll(k2)
integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1)
c0 = SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)
c1 = SCF_density_matrix_ao_alpha(k,i)
c2 = SCF_density_matrix_ao_beta(k,i)
if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then
cycle
endif
if (values(1) == 0.d0) then
values(1) = ao_two_e_integral(k0,l0,i0,j0)
endif
integral = c0 * values(1)
ao_two_e_integral_alpha_tmp(i,j) += integral
ao_two_e_integral_beta_tmp (i,j) += integral
integral = values(k1)
ao_two_e_integral_alpha_tmp(l,j) -= SCF_density_matrix_ao_alpha(k,i) * integral
ao_two_e_integral_beta_tmp (l,j) -= SCF_density_matrix_ao_beta (k,i) * integral
integral = values(1)
ao_two_e_integral_alpha_tmp(l,j) -= c1 * integral
ao_two_e_integral_beta_tmp (l,j) -= c2 * integral
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp
ao_two_e_integral_beta += ao_two_e_integral_beta_tmp
!$OMP END CRITICAL
deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)
!$OMP END PARALLEL
!$OMP END DO NOWAIT
!$OMP CRITICAL
ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp
ao_two_e_integral_beta += ao_two_e_integral_beta_tmp
!$OMP END CRITICAL
deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)
!$OMP END PARALLEL
else
PROVIDE ao_two_e_integrals_in_map
integer(omp_lock_kind) :: lck(ao_num)
integer(map_size_kind) :: i8
integer :: ii(8), jj(8), kk(8), ll(8), k2
integer(cache_map_size_kind) :: n_elements_max, n_elements
integer(key_kind), allocatable :: keys(:)
double precision, allocatable :: values(:)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max,&
!$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)&
!$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
!$OMP ao_integrals_map, ao_two_e_integral_alpha, ao_two_e_integral_beta)
call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max)
allocate(keys(n_elements_max), values(n_elements_max))
allocate(ao_two_e_integral_alpha_tmp(ao_num,ao_num), &
ao_two_e_integral_beta_tmp(ao_num,ao_num))
ao_two_e_integral_alpha_tmp = 0.d0
ao_two_e_integral_beta_tmp = 0.d0
!$OMP DO SCHEDULE(static,1)
do i8=0_8,ao_integrals_map%map_size
n_elements = n_elements_max
call get_cache_map(ao_integrals_map,i8,keys,values,n_elements)
do k1=1,n_elements
call two_e_integrals_index_reverse(kk,ii,ll,jj,keys(k1))
do k2=1,8
if (kk(k2)==0) then
cycle
endif
i = ii(k2)
j = jj(k2)
k = kk(k2)
l = ll(k2)
integral = (SCF_density_matrix_ao_alpha(k,l)+SCF_density_matrix_ao_beta(k,l)) * values(k1)
ao_two_e_integral_alpha_tmp(i,j) += integral
ao_two_e_integral_beta_tmp (i,j) += integral
integral = values(k1)
ao_two_e_integral_alpha_tmp(l,j) -= SCF_density_matrix_ao_alpha(k,i) * integral
ao_two_e_integral_beta_tmp (l,j) -= SCF_density_matrix_ao_beta (k,i) * integral
enddo
enddo
enddo
!$OMP END DO NOWAIT
!$OMP CRITICAL
ao_two_e_integral_alpha += ao_two_e_integral_alpha_tmp
ao_two_e_integral_beta += ao_two_e_integral_beta_tmp
!$OMP END CRITICAL
deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)
!$OMP END PARALLEL
endif
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_two_e_integral_alpha_chol, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta_chol , (ao_num, ao_num) ]
use map_module
implicit none
BEGIN_DOC
! Alpha and Beta Fock matrices in AO basis set
END_DOC
integer :: m,n,l,s,j
double precision :: integral
double precision, allocatable :: X(:), X2(:,:,:,:), X3(:,:,:,:)
allocate (X(cholesky_ao_num))
! X(j) = \sum_{mn} SCF_density_matrix_ao(m,n) * cholesky_ao(m,n,j)
call dgemm('T','N',cholesky_ao_num,1,ao_num*ao_num,1.d0, &
cholesky_ao, ao_num*ao_num, &
SCF_density_matrix_ao, ao_num*ao_num,0.d0, &
X, cholesky_ao_num)
!
! ao_two_e_integral_alpha(m,n) = \sum_{j} cholesky_ao(m,n,j) * X(j)
call dgemm('N','N',ao_num*ao_num,1,cholesky_ao_num, 1.d0, &
cholesky_ao, ao_num*ao_num, &
X, cholesky_ao_num, 0.d0, &
ao_two_e_integral_alpha_chol, ao_num*ao_num)
deallocate(X)
ao_two_e_integral_beta_chol = ao_two_e_integral_alpha_chol
allocate(X2(ao_num,ao_num,cholesky_ao_num,2))
! ao_two_e_integral_alpha_chol (l,s) -= cholesky_ao(l,m,j) * SCF_density_matrix_ao_beta (m,n) * cholesky_ao(n,s,j)
call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0, &
SCF_density_matrix_ao_alpha, ao_num, &
cholesky_ao, ao_num, 0.d0, &
X2(1,1,1,1), ao_num)
call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0, &
SCF_density_matrix_ao_beta, ao_num, &
cholesky_ao, ao_num, 0.d0, &
X2(1,1,1,2), ao_num)
allocate(X3(ao_num,cholesky_ao_num,ao_num,2))
do s=1,ao_num
do j=1,cholesky_ao_num
do m=1,ao_num
X3(m,j,s,1) = X2(m,s,j,1)
X3(m,j,s,2) = X2(m,s,j,2)
enddo
enddo
enddo
deallocate(X2)
call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0, &
cholesky_ao, ao_num, &
X3(1,1,1,1), ao_num*cholesky_ao_num, 1.d0, &
ao_two_e_integral_alpha_chol, ao_num)
call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0, &
cholesky_ao, ao_num, &
X3(1,1,1,2), ao_num*cholesky_ao_num, 1.d0, &
ao_two_e_integral_beta_chol, ao_num)
deallocate(X3)
END_PROVIDER
BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ]

9
src/hartree_fock/scf.irp.f
View File

@ -80,9 +80,14 @@ subroutine run
mo_label = 'Orthonormalized'
call Roothaan_Hall_SCF
call ezfio_set_hartree_fock_energy(SCF_energy)
write(json_unit,json_array_open_fmt) 'scf'
call Roothaan_Hall_SCF
write(json_unit,json_array_close_fmtx)
call json_close
call ezfio_set_hartree_fock_energy(SCF_energy)
end

24
src/iterations/EZFIO.cfg
View File

@ -1,24 +0,0 @@
[n_iter]
interface: ezfio
doc: Number of saved iterations
type:integer
default: 1
[n_det_iterations]
interface: ezfio, provider
doc: Number of determinants at each iteration
type: integer
size: (100)
[energy_iterations]
interface: ezfio, provider
doc: The variational energy at each iteration
type: double precision
size: (determinants.n_states,100)
[pt2_iterations]
interface: ezfio, provider
doc: The |PT2| correction at each iteration
type: double precision
size: (determinants.n_states,100)

37
src/iterations/io.irp.f
View File

@ -1,37 +0,0 @@
BEGIN_PROVIDER [ integer, n_iter ]
implicit none
BEGIN_DOC
! number of iterations
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
double precision :: zeros(N_states,100)
integer :: izeros(100)
zeros = 0.d0
izeros = 0
call ezfio_set_iterations_n_iter(0)
call ezfio_set_iterations_energy_iterations(zeros)
call ezfio_set_iterations_pt2_iterations(zeros)
call ezfio_set_iterations_n_det_iterations(izeros)
n_iter = 1
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( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read n_iter with MPI'
endif
IRP_ENDIF
call write_time(6)
END_PROVIDER

91
src/iterations/iterations.irp.f
View File

@ -1,42 +1,65 @@
BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ]
implicit none
BEGIN_DOC
! Extrapolated energy, using E_var = f(PT2) where PT2=0
END_DOC
integer :: i
do i=1,min(N_states,N_det)
call extrapolate_data(N_iter, &
energy_iterations(i,1:N_iter), &
pt2_iterations(i,1:N_iter), &
extrapolated_energy(1:N_iter,i))
enddo
END_PROVIDER
subroutine save_iterations(e_, pt2_,n_)
BEGIN_PROVIDER [ integer, N_iter ]
implicit none
BEGIN_DOC
! Update the energy in the EZFIO file.
! Number of CIPSI iterations
END_DOC
integer, intent(in) :: n_
double precision, intent(in) :: e_(N_states), pt2_(N_states)
integer :: i
if (N_iter == 101) then
do i=2,N_iter-1
energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter)
pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter)
N_iter = 0
END_PROVIDER
BEGIN_PROVIDER [ integer, N_iter_max ]
implicit none
BEGIN_DOC
! Max number of iterations for extrapolations
END_DOC
N_iter_max = 8
END_PROVIDER
BEGIN_PROVIDER [ double precision, energy_iterations , (n_states,N_iter_max) ]
&BEGIN_PROVIDER [ double precision, pt2_iterations , (n_states,N_iter_max) ]
&BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter_max,N_states) ]
implicit none
BEGIN_DOC
! The energy at each iteration for the extrapolations
END_DOC
energy_iterations = 0.d0
pt2_iterations = 0.d0
extrapolated_energy = 0.d0
END_PROVIDER
subroutine increment_n_iter(e, pt2_data)
use selection_types
implicit none
BEGIN_DOC
! Does what is necessary to increment n_iter
END_DOC
double precision, intent(in) :: e(*)
type(pt2_type), intent(in) :: pt2_data
integer :: k, i
if (N_det < N_states) return
if (N_iter < N_iter_max) then
N_iter += 1
else
do k=2,N_iter
energy_iterations(1:N_states,k-1) = energy_iterations(1:N_states,k)
pt2_iterations(1:N_states,k-1) = pt2_iterations(1:N_states,k)
enddo
N_iter = N_iter-1
TOUCH N_iter
endif
energy_iterations(1:N_states,N_iter) = e(1:N_states)
pt2_iterations(1:N_states,N_iter) = pt2_data % rpt2(1:N_states)
energy_iterations(1:N_states,N_iter) = e_(1:N_states)
pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states)
n_det_iterations(N_iter) = n_
call ezfio_set_iterations_N_iter(N_iter)
call ezfio_set_iterations_energy_iterations(energy_iterations)
call ezfio_set_iterations_pt2_iterations(pt2_iterations)
call ezfio_set_iterations_n_det_iterations(n_det_iterations)
if (N_iter < 2) then
extrapolated_energy(1,:) = energy_iterations(:,1) + pt2_iterations(:,1)
extrapolated_energy(2,:) = energy_iterations(:,2) + pt2_iterations(:,2)
else
do i=1,N_states
call extrapolate_data(N_iter, &
energy_iterations(i,1:N_iter), &
pt2_iterations(i,1:N_iter), &
extrapolated_energy(1:N_iter,i))
enddo
endif
end

14
src/iterations/print_extrapolation.irp.f
View File

@ -5,10 +5,14 @@ subroutine print_extrapolated_energy
END_DOC
integer :: i,k
integer :: N_states_p, N_iter_p
if (N_iter< 2) then
return
endif
N_states_p = min(N_states,N_det)
N_iter_p = min(N_iter, 8)
write(*,'(A)') ''
write(*,'(A)') 'Extrapolated energies'
write(*,'(A)') '------------------------'
@ -20,20 +24,20 @@ subroutine print_extrapolated_energy
write(*,*) '=========== ', '==================='
write(*,*) 'minimum PT2 ', 'Extrapolated energy'
write(*,*) '=========== ', '==================='
do k=2,min(N_iter,8)
write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1)
do k=2,N_iter_p
write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,k), extrapolated_energy(k,1)
enddo
write(*,*) '=========== ', '==================='
do i=2, min(N_states,N_det)
do i=2, N_states_p
print *, ''
print *, 'State ', i
print *, ''
write(*,*) '=========== ', '=================== ', '=================== ', '==================='
write(*,*) 'minimum PT2 ', 'Extrapolated energy ', ' Excitation (a.u) ', ' Excitation (eV) '
write(*,*) '=========== ', '=================== ', '=================== ', '==================='
do k=2,min(N_iter,8)
write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
do k=2,N_iter_p
write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,k), extrapolated_energy(k,i), &
extrapolated_energy(k,i) - extrapolated_energy(k,1), &
(extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
enddo

24
src/iterations_tc/EZFIO.cfg
View File

@ -1,24 +0,0 @@
[n_iter]
interface: ezfio
doc: Number of saved iterations
type:integer
default: 1
[n_det_iterations]
interface: ezfio, provider
doc: Number of determinants at each iteration
type: integer
size: (100)
[energy_iterations]
interface: ezfio, provider
doc: The variational energy at each iteration
type: double precision
size: (determinants.n_states,100)
[pt2_iterations]
interface: ezfio, provider
doc: The |PT2| correction at each iteration
type: double precision
size: (determinants.n_states,100)

0
src/iterations_tc/NEED
View File

37
src/iterations_tc/io.irp.f
View File

@ -1,37 +0,0 @@
BEGIN_PROVIDER [ integer, n_iter ]
implicit none
BEGIN_DOC
! number of iterations
END_DOC
logical :: has
PROVIDE ezfio_filename
if (mpi_master) then
double precision :: zeros(N_states,100)
integer :: izeros(100)
zeros = 0.d0
izeros = 0
call ezfio_set_iterations_n_iter(0)
call ezfio_set_iterations_energy_iterations(zeros)
call ezfio_set_iterations_pt2_iterations(zeros)
call ezfio_set_iterations_n_det_iterations(izeros)
n_iter = 1
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( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read n_iter with MPI'
endif
IRP_ENDIF
call write_time(6)
END_PROVIDER

43
src/iterations_tc/iterations.irp.f
View File

@ -1,43 +0,0 @@
BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ]
implicit none
BEGIN_DOC
! Extrapolated energy, using E_var = f(PT2) where PT2=0
END_DOC
! integer :: i
extrapolated_energy = 0.D0
END_PROVIDER
subroutine get_extrapolated_energy(Niter,ept2,pt1,extrap_energy)
implicit none
integer, intent(in) :: Niter
double precision, intent(in) :: ept2(Niter),pt1(Niter),extrap_energy(Niter)
call extrapolate_data(Niter,ept2,pt1,extrap_energy)
end
subroutine save_iterations(e_, pt2_,n_)
implicit none
BEGIN_DOC
! Update the energy in the EZFIO file.
END_DOC
integer, intent(in) :: n_
double precision, intent(in) :: e_(N_states), pt2_(N_states)
integer :: i
if (N_iter == 101) then
do i=2,N_iter-1
energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter)
pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter)
enddo
N_iter = N_iter-1
TOUCH N_iter
endif
energy_iterations(1:N_states,N_iter) = e_(1:N_states)
pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states)
n_det_iterations(N_iter) = n_
call ezfio_set_iterations_N_iter(N_iter)
call ezfio_set_iterations_energy_iterations(energy_iterations)
call ezfio_set_iterations_pt2_iterations(pt2_iterations)
call ezfio_set_iterations_n_det_iterations(n_det_iterations)
end

46
src/iterations_tc/print_extrapolation.irp.f
View File

@ -1,46 +0,0 @@
subroutine print_extrapolated_energy
implicit none
BEGIN_DOC
! Print the extrapolated energy in the output
END_DOC
integer :: i,k
if (N_iter< 2) then
return
endif
write(*,'(A)') ''
write(*,'(A)') 'Extrapolated energies'
write(*,'(A)') '------------------------'
write(*,'(A)') ''
print *, ''
print *, 'State ', 1
print *, ''
write(*,*) '=========== ', '==================='
write(*,*) 'minimum PT2 ', 'Extrapolated energy'
write(*,*) '=========== ', '==================='
do k=2,min(N_iter,8)
write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1)
enddo
write(*,*) '=========== ', '==================='
do i=2, min(N_states,N_det)
print *, ''
print *, 'State ', i
print *, ''
write(*,*) '=========== ', '=================== ', '=================== ', '==================='
write(*,*) 'minimum PT2 ', 'Extrapolated energy ', ' Excitation (a.u) ', ' Excitation (eV) '
write(*,*) '=========== ', '=================== ', '=================== ', '==================='
do k=2,min(N_iter,8)
write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
extrapolated_energy(k,i) - extrapolated_energy(k,1), &
(extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
enddo
write(*,*) '=========== ', '=================== ', '=================== ', '==================='
enddo
print *, ''
end subroutine

104
src/iterations_tc/print_summary.irp.f
View File

@ -1,104 +0,0 @@
subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s2_)
use selection_types
implicit none
BEGIN_DOC
! Print the extrapolated energy in the output
END_DOC
integer, intent(in) :: n_det_, n_configuration_, n_st
double precision, intent(in) :: e_(n_st), s2_(n_st)
type(pt2_type) , intent(in) :: pt2_data, pt2_data_err
integer :: i, k
integer :: N_states_p
character*(9) :: pt2_string
character*(512) :: fmt
if (do_pt2) then
pt2_string = ' '
else
pt2_string = '(approx)'
endif
N_states_p = min(N_det_,n_st)
print *, ''
print '(A,I12)', 'Summary at N_det = ', N_det_
print '(A)', '-----------------------------------'
print *, ''
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))'
write(*,fmt) ('State',k, k=1,N_states_p)
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))'
write(*,fmt) '# E ', e_(1:N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
endif
write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p)
write(*,fmt) '# rPT2'//pt2_string, (pt2_data % rpt2(k), pt2_data_err % rpt2(k), k=1,N_states_p)
write(*,'(A)') '#'
write(*,fmt) '# E+PT2 ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p)
write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_data % rpt2(k),pt2_data_err % rpt2(k), k=1,N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), &
dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p)
write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, &
dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p)
endif
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
print *, ''
print *, 'N_det = ', N_det_
print *, 'N_states = ', n_st
if (s2_eig) then
print *, 'N_cfg = ', N_configuration_
if (only_expected_s2) then
print *, 'N_csf = ', N_csf
endif
endif
print *, ''
do k=1, N_states_p
print*,'* State ',k
print *, '< S^2 > = ', s2_(k)
print *, 'E = ', e_(k)
print *, 'Variance = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k)
print *, 'PT norm = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k))
print *, 'PT2 = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
print *, 'rPT2 = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k)
print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k)
print *, ''
enddo
print *, '-----'
if(n_st.gt.1)then
print *, 'Variational Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (e_(i) - e_(1)), &
(e_(i) - e_(1)) * 27.211396641308d0
enddo
print *, '-----'
print*, 'Variational + perturbative Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), &
(e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0
enddo
print *, '-----'
print*, 'Variational + renormalized perturbative Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), &
(e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0
enddo
endif
! call print_energy_components()
end subroutine

5
src/json/EZFIO.cfg Normal file
View File

@ -0,0 +1,5 @@
[empty]
type: logical
doc: Needed to create the json directory
interface: ezfio

1
src/json/NEED Normal file
View File

@ -0,0 +1 @@
ezfio_files

5
src/json/README.rst Normal file
View File

@ -0,0 +1,5 @@
====
json
====
JSON files to simplify getting output information from QP.

45
src/json/json.irp.f Normal file
View File

@ -0,0 +1,45 @@
BEGIN_PROVIDER [ character*(128), json_filename ]
implicit none
BEGIN_DOC
! Fortran unit of the JSON file
END_DOC
integer, external :: getUnitAndOpen
integer :: counter
character*(128) :: prefix
logical :: exists
prefix = trim(ezfio_filename)//'/json/'
call lock_io
exists = .True.
counter = 0
do while (exists)
counter += 1
write(json_filename, '(A,I5.5,A)') trim(prefix), counter, '.json'
INQUIRE(FILE=trim(json_filename), EXIST=exists)
enddo
call unlock_io
END_PROVIDER
BEGIN_PROVIDER [ integer, json_unit]
implicit none
BEGIN_DOC
! Unit file for JSON output
END_DOC
integer, external :: getUnitAndOpen
call ezfio_set_json_empty(.False.)
call lock_io
json_unit = getUnitAndOpen(json_filename, 'w')
write(json_unit, '(A)') '{'
call unlock_io
END_PROVIDER
subroutine json_close
call lock_io
write(json_unit, '(A)') '}'
close(json_unit)
call unlock_io
FREE json_unit
end

46
src/json/json_formats.irp.f Normal file
View File

@ -0,0 +1,46 @@
BEGIN_PROVIDER [ character*(64), json_int_fmt ]
&BEGIN_PROVIDER [ character*(64), json_int_fmtx ]
&BEGIN_PROVIDER [ character*(64), json_real_fmt ]
&BEGIN_PROVIDER [ character*(64), json_real_fmtx ]
&BEGIN_PROVIDER [ character*(64), json_str_fmt ]
&BEGIN_PROVIDER [ character*(64), json_str_fmtx ]
&BEGIN_PROVIDER [ character*(64), json_true_fmt ]
&BEGIN_PROVIDER [ character*(64), json_true_fmtx ]
&BEGIN_PROVIDER [ character*(64), json_false_fmt ]
&BEGIN_PROVIDER [ character*(64), json_false_fmtx ]
&BEGIN_PROVIDER [ character*(64), json_array_open_fmt ]
&BEGIN_PROVIDER [ character*(64), json_array_uopen_fmt ]
&BEGIN_PROVIDER [ character*(64), json_array_close_fmt ]
&BEGIN_PROVIDER [ character*(64), json_array_close_uopen_fmt ]
&BEGIN_PROVIDER [ character*(64), json_array_close_fmtx ]
&BEGIN_PROVIDER [ character*(64), json_dict_open_fmt ]
&BEGIN_PROVIDER [ character*(64), json_dict_uopen_fmt ]
&BEGIN_PROVIDER [ character*(64), json_dict_close_uopen_fmt ]
&BEGIN_PROVIDER [ character*(64), json_dict_close_fmt ]
&BEGIN_PROVIDER [ character*(64), json_dict_close_fmtx ]
implicit none
BEGIN_DOC
! Formats for JSON output.
! x: used to mark the last write (no comma)
END_DOC
json_int_fmt = '('' "'',A,''": '',I10,'','')'
json_int_fmtx = '('' "'',A,''": '',I10)'
json_real_fmt = '('' "'',A,''": '',E22.15,'','')'
json_real_fmtx = '('' "'',A,''": '',E22.15)'
json_str_fmt = '('' "'',A,''": "'',A,''",'')'
json_str_fmtx = '('' "'',A,''": "'',A,''"'')'
json_true_fmt = '('' "'',A,''": true,'')'
json_true_fmtx = '('' "'',A,''": true'')'
json_false_fmt = '('' "'',A,''": false,'')'
json_false_fmtx = '('' "'',A,''": false'')'
json_array_open_fmt = '('' "'',A,''": ['')'
json_array_uopen_fmt = '('' ['')'
json_array_close_fmt = '('' ],'')'
json_array_close_uopen_fmt = '('' ], ['')'
json_array_close_fmtx = '('' ]'')'
json_dict_open_fmt = '('' "'',A,''": {'')'
json_dict_uopen_fmt = '('' {'')'
json_dict_close_fmt = '('' },'')'
json_dict_close_uopen_fmt = '('' }, {'')'
json_dict_close_fmtx = '('' }'')'
END_PROVIDER

4
src/kohn_sham/ks_scf.irp.f
View File

@ -90,7 +90,11 @@ subroutine run
! Choose SCF algorithm
write(json_unit,*) '"scf" : ['
call Roothaan_Hall_SCF
write(json_unit,*) ']'
call json_close
end

3
src/kohn_sham_rs/rs_ks_scf.irp.f
View File

@ -93,7 +93,10 @@ subroutine run
level_shift += 1.d0
touch level_shift
write(json_unit,*) '"scf" : ['
call Roothaan_Hall_SCF
write(json_unit,*) ']'
call json_close
call ezfio_set_kohn_sham_rs_energy(SCF_energy)
write(*, '(A22,X,F16.10)') 'one_e_energy = ',one_e_energy

16
src/mo_two_e_ints/cholesky.irp.f Normal file
View File

@ -0,0 +1,16 @@
BEGIN_PROVIDER [ double precision, cholesky_mo, (mo_num, mo_num, cholesky_ao_num) ]
implicit none
BEGIN_DOC
! Cholesky vectors in MO basis
END_DOC
integer :: k
!$OMP PARALLEL DO PRIVATE(k)
do k=1,cholesky_ao_num
call ao_to_mo(cholesky_ao(1,1,k),ao_num,cholesky_mo(1,1,k),mo_num)
enddo
!$OMP END PARALLEL DO
END_PROVIDER

2
src/mo_two_e_ints/map_integrals.irp.f
View File

@ -377,6 +377,7 @@ integer function load_mo_integrals(filename)
integer*8 :: n, j
load_mo_integrals = 1
open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
call lock_io()
read(66,err=98,end=98) iknd, kknd
if (iknd /= integral_kind) then
print *, 'Wrong integrals kind in file :', iknd
@ -399,6 +400,7 @@ integer function load_mo_integrals(filename)
n = mo_integrals_map%map(i)%n_elements
read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
enddo
call unlock_io()
call map_sort(mo_integrals_map)
load_mo_integrals = 0
return

74
src/mo_two_e_ints/mo_bi_integrals.irp.f
View File

@ -50,13 +50,16 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
call cpu_time(cpu_1)
if(no_vvvv_integrals)then
! call four_idx_novvvv
call four_idx_novvvv_old
else
if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then
call four_idx_dgemm
if (do_ao_cholesky) then
call add_integrals_to_map_cholesky
else
call add_integrals_to_map(full_ijkl_bitmask_4)
if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then
call four_idx_dgemm
else
call add_integrals_to_map(full_ijkl_bitmask_4)
endif
endif
endif
@ -175,7 +178,7 @@ subroutine add_integrals_to_map(mask_ijkl)
implicit none
BEGIN_DOC
! Adds integrals to tha MO map according to some bitmask
! Adds integrals to the MO map according to some bitmask
END_DOC
integer(bit_kind), intent(in) :: mask_ijkl(N_int,4)
@ -450,13 +453,72 @@ subroutine add_integrals_to_map(mask_ijkl)
end
subroutine add_integrals_to_map_cholesky
use bitmasks
implicit none
BEGIN_DOC
! Adds integrals to the MO map using Cholesky vectors
END_DOC
integer :: i,j,k,l,m
integer :: size_buffer, n_integrals
size_buffer = min(mo_num*mo_num*mo_num,16000000)
double precision, allocatable :: Vtmp(:,:,:,:)
integer(key_kind) , allocatable :: buffer_i(:)
real(integral_kind), allocatable :: buffer_value(:)
if (.True.) then
! In-memory transformation
allocate (Vtmp(mo_num,mo_num,mo_num,mo_num))
call dgemm('N','T',mo_num*mo_num,mo_num*mo_num,cholesky_ao_num,1.d0, &
cholesky_mo, mo_num*mo_num, &
cholesky_mo, mo_num*mo_num, 0.d0, &
Vtmp, mo_num*mo_num)
!$OMP PARALLEL PRIVATE(i,j,k,l,n_integrals,buffer_value, buffer_i)
allocate (buffer_i(size_buffer), buffer_value(size_buffer))
n_integrals = 0
!$OMP DO
do l=1,mo_num
do k=1,l
do j=1,mo_num
do i=1,j
if (abs(Vtmp(i,j,k,l)) > mo_integrals_threshold) then
n_integrals += 1
buffer_value(n_integrals) = Vtmp(i,j,k,l)
!DIR$ FORCEINLINE
call mo_two_e_integrals_index(i,k,j,l,buffer_i(n_integrals))
if (n_integrals == size_buffer) then
call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
n_integrals = 0
endif
endif
enddo
enddo
enddo
enddo
!$OMP END DO
call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
deallocate(buffer_i, buffer_value)
!$OMP END PARALLEL
deallocate(Vtmp)
call map_unique(mo_integrals_map)
endif
end
subroutine add_integrals_to_map_three_indices(mask_ijk)
use bitmasks
implicit none
BEGIN_DOC
! Adds integrals to tha MO map according to some bitmask
! Adds integrals to the MO map according to some bitmask
END_DOC
integer(bit_kind), intent(in) :: mask_ijk(N_int,3)

14
src/nuclei/nuclei.irp.f
View File

@ -241,13 +241,13 @@ END_PROVIDER
enddo
character*(80) :: buffer, dummy
do
read(iunit,'(A80)',end=10) buffer
read(buffer,*) i ! First read i
read(buffer,*) i, element_name(i), dummy, element_mass(i)
enddo
10 continue
close(10)
endif
read(iunit,'(A80)',end=10) buffer
read(buffer,*) i ! First read i
read(buffer,*) i, element_name(i), dummy, element_mass(i)
enddo
10 continue
close(10)
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank

1
src/scf_utils/NEED
View File

@ -1,2 +1,3 @@
mo_guess
bitmask
json

43
src/scf_utils/roothaan_hall_scf.irp.f
View File

@ -12,6 +12,7 @@ END_DOC
integer :: iteration_SCF,dim_DIIS,index_dim_DIIS
logical :: converged
integer :: i,j
logical, external :: qp_stop
double precision, allocatable :: mo_coef_save(:,:)
@ -50,10 +51,8 @@ END_DOC
!
PROVIDE FPS_SPF_matrix_AO Fock_matrix_AO
do while ( &
( (max_error_DIIS > threshold_DIIS_nonzero) .or. &
(dabs(Delta_energy_SCF) > thresh_SCF) &
) .and. (iteration_SCF < n_it_SCF_max) )
converged = .False.
do while ( .not.converged .and. (iteration_SCF < n_it_SCF_max) )
! Increment cycle number
@ -144,17 +143,49 @@ END_DOC
SOFT_TOUCH level_shift
energy_SCF_previous = energy_SCF
converged = ( (max_error_DIIS <= threshold_DIIS_nonzero) .and. &
(dabs(Delta_energy_SCF) <= thresh_SCF) )
! Print results at the end of each iteration
write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') &
iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS
! Write data in JSON file
call lock_io
if (iteration_SCF == 1) then
write(json_unit, json_dict_uopen_fmt)
else
write(json_unit, json_dict_close_uopen_fmt)
endif
write(json_unit, json_int_fmt) 'iteration', iteration_SCF
write(json_unit, json_real_fmt) 'energy', energy_SCF
write(json_unit, json_real_fmt) 'delta_energy_SCF', Delta_energy_SCF
write(json_unit, json_real_fmt) 'max_error_DIIS', max_error_DIIS
write(json_unit, json_real_fmt) 'level_shift', level_shift
write(json_unit, json_int_fmt) 'dim_DIIS', dim_DIIS
call unlock_io
if (Delta_energy_SCF < 0.d0) then
call save_mos
write(json_unit, json_true_fmt) 'saved'
else
write(json_unit, json_false_fmt) 'saved'
endif
call lock_io
if (converged) then
write(json_unit, json_true_fmtx) 'converged'
else
write(json_unit, json_false_fmtx) 'converged'
endif
call unlock_io
if (qp_stop()) exit
enddo
write(json_unit, json_dict_close_fmtx)
if (iteration_SCF < n_it_SCF_max) then
mo_label = 'Canonical'
@ -166,6 +197,10 @@ END_DOC
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
'====','================','================','================','================'
write(6,*)
if (converged) then
write(6,*) 'SCF converged'
endif
if(.not.frozen_orb_scf)then
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1), &

10
src/tc_bi_ortho/h_tc_s2_u0.irp.f
View File

@ -152,7 +152,7 @@ subroutine H_tc_s2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishi
use bitmasks
implicit none
BEGIN_DOC
! Computes $v_t = H | u_t\rangle$
! Computes $v_t = H | u_t\rangle$
!
! Default should be 1,N_det,0,1
!
@ -160,7 +160,7 @@ subroutine H_tc_s2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishi
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
logical, intent(in) :: do_right
logical, intent(in) :: do_right
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
@ -193,7 +193,7 @@ subroutine H_tc_s2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,ie
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
logical, intent(in) :: do_right
logical, intent(in) :: do_right
double precision, intent(out) :: v_t(N_st,sze), s_t(N_st,sze)
double precision :: hij, sij
@ -570,7 +570,7 @@ compute_singles=.True.
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
tmp_det2(1:N_int,1) = psi_det_alpha_unique(1:N_int, lrow)
! call i_H_j( tmp_det, tmp_det2, $N_int, hij)
! call i_H_j_double_spin( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij)
if(do_right)then
@ -721,7 +721,7 @@ compute_singles=.True.
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol)
tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol)
! call i_H_j( tmp_det, tmp_det2, $N_int, hij)
! call i_H_j_double_spin( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij)
if(do_right)then

80
src/tools/print_e_conv.irp.f
View File

@ -1,80 +0,0 @@
program print_e_conv
implicit none
BEGIN_DOC
! program that prints in a human readable format the convergence of the CIPSI algorithm.
!
! for all istate, this program produces
!
! * a file "EZFIO.istate.conv" containing the variational and var+PT2 energies as a function of N_det
!
! * for istate > 1, a file EZFIO.istate.delta_e.conv containing the energy difference (both var and var+PT2) with the ground state as a function of N_det
END_DOC
provide ezfio_filename
call routine_e_conv
end
subroutine routine_e_conv
implicit none
BEGIN_DOC
! routine called by :c:func:`print_e_conv`
END_DOC
integer :: N_iter_tmp
integer :: i,istate
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
character*(128) :: filename
integer, allocatable :: n_det_tmp(:)
call ezfio_get_iterations_N_iter(N_iter_tmp)
print*,'N_iter_tmp = ',N_iter_tmp
double precision, allocatable :: e(:,:),pt2(:,:)
allocate(e(N_states, 100),pt2(N_states, 100),n_det_tmp(100))
call ezfio_get_iterations_energy_iterations(e)
call ezfio_get_iterations_pt2_iterations(pt2)
call ezfio_get_iterations_n_det_iterations(n_det_tmp)
do istate = 1, N_states
if (istate.lt.10)then
write (filename, "(I1)")istate
else
write (filename, "(I2)")istate
endif
print*,filename
output=trim(ezfio_filename)//'.'//trim(filename)//'.conv'
output=trim(output)
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,*)'# N_det E_var E_var + PT2'
do i = 1, N_iter_tmp
write(i_unit_output,'(I9,X,3(F16.10,X))')n_det_tmp(i),e(istate,i),e(istate,i) + pt2(istate,i)
enddo
enddo
if(N_states.gt.1)then
double precision, allocatable :: deltae(:,:),deltae_pt2(:,:)
allocate(deltae(N_states,100),deltae_pt2(N_states,100))
do i = 1, N_iter_tmp
do istate = 1, N_states
deltae(istate,i) = e(istate,i) - e(1,i)
deltae_pt2(istate,i) = e(istate,i) + pt2(istate,i) - (e(1,i) + pt2(1,i))
enddo
enddo
do istate = 2, N_states
if (istate.lt.10)then
write (filename, "(I1)")istate
else
write (filename, "(I2)")istate
endif
output=trim(ezfio_filename)//'.'//trim(filename)//'.delta_e.conv'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,*)'# N_det Delta E_var Delta (E_var + PT2)'
do i = 1, N_iter_tmp
write(i_unit_output,'(I9,X,100(F16.10,X))')n_det_tmp(i),deltae(istate,i),deltae_pt2(istate,i)
enddo
enddo
endif
end

12
src/two_body_rdm/io_two_rdm.irp.f
View File

@ -1,15 +1,17 @@
subroutine write_array_two_rdm(n_orb,nstates,array_tmp,name_file)
implicit none
integer, intent(in) :: n_orb,nstates
character*(128), intent(in) :: name_file
character*(128), intent(in) :: name_file
double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'W')
call lock_io()
write(i_unit_output)array_tmp
call unlock_io()
close(unit=i_unit_output)
end
@ -18,12 +20,14 @@ subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer, intent(in) :: n_orb,nstates
character*(128), intent(in) :: name_file
character*(128), intent(in) :: name_file
double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states)
PROVIDE ezfio_filename
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'R')
call lock_io()
read(i_unit_output)array_tmp
call unlock_io()
close(unit=i_unit_output)
end

16
src/utils/format_w_error.irp.f
View File

@ -1,7 +1,7 @@
subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error)
implicit none
BEGIN_DOC
! Format for double precision, value(error)
END_DOC
@ -14,7 +14,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err
! out
! | format_value | character | string FX.Y for the format |
! | str_error | character | string of the error |
! | str_error | character | string of the error |
! internal
! | str_size | character | size in string format |
@ -33,6 +33,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err
character(len=20) :: str_size, str_nb_digits, str_exp
integer :: nb_digits
call lock_io()
! max_nb_digit: Y max
! size_nb = Size of the double: X (FX.Y)
write(str_size,'(I3)') size_nb
@ -40,17 +41,17 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err
! Error
write(str_exp,'(1pE20.0)') error
str_error = trim(adjustl(str_exp))
! Number of digit: Y (FX.Y) from the exponent
str_nb_digits = str_exp(19:20)
read(str_nb_digits,*) nb_digits
! If the error is 0d0
if (error <= 1d-16) then
if (error <= 1d-16) then
write(str_nb_digits,*) max_nb_digits
endif
! If the error is too small
! If the error is too small
if (nb_digits > max_nb_digits) then
write(str_nb_digits,*) max_nb_digits
str_error(1:1) = '0'
@ -65,7 +66,8 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err
! FX.Y,A1,A1,A1 for value(str_error)
!string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1'
! FX.Y just for the value
! FX.Y just for the value
format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))
call unlock_io()
end

2
src/utils/linear_algebra.irp.f
View File

@ -1854,7 +1854,7 @@ do k = 1, N
end do
! TODO: It should be possible to use only one vector of size (1:rank) as a buffer
! to do the swapping in-place
U = 0.00D+0
U(:,:) = 0.00D+0
do k = 1, N
l = piv(k)
U(l, :) = A(1:rank, k)

8
src/utils/memory.irp.f
View File

@ -8,7 +8,9 @@ BEGIN_PROVIDER [ integer, qp_max_mem ]
qp_max_mem = 2000
call getenv('QP_MAXMEM',env)
if (trim(env) /= '') then
call lock_io()
read(env,*) qp_max_mem
call unlock_io()
endif
call write_int(6,qp_max_mem,'Target maximum memory (GB)')
@ -25,7 +27,7 @@ subroutine resident_memory(value)
character*(32) :: key
double precision, intent(out) :: value
call omp_set_lock(file_lock)
call lock_io()
call usleep(10)
value = 0.d0
@ -40,7 +42,7 @@ subroutine resident_memory(value)
20 continue
close(iunit)
value = value / (1024.d0*1024.d0)
call omp_unset_lock(file_lock)
call unlock_io()
end function
subroutine total_memory(value)
@ -53,6 +55,7 @@ subroutine total_memory(value)
character*(32) :: key
double precision, intent(out) :: value
call lock_io()
iunit = getUnitAndOpen('/proc/self/status','r')
do
read(iunit,*,err=10,end=20) key, value
@ -64,6 +67,7 @@ subroutine total_memory(value)
20 continue
close(iunit)
value = value / (1024.d0*1024.d0)
call unlock_io()
end function
double precision function memory_of_double(n)