diff --git a/README.md b/README.md
index b8141b88..24d6277e 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,7 @@
# Quantum Package 2.2
-
+
+
[](https://zenodo.org/badge/latestdoi/167513335)
@@ -44,9 +45,19 @@ https://arxiv.org/abs/1902.08154
# Credits
+* [TREX Center of Excellence](https://trex-coe.eu)
+* [ERC PTEROSOR](https://lcpq.github.io/PTEROSOR)
* [CNRS](http://www.cnrs.fr)
* [Laboratoire de Chimie et Physique Quantiques](http://lcpq.ups-tlse.fr)
* [Laboratoire de Chimie Théorique](http://www.lct.jussieu.fr)
* [Argonne Leadership Computing Facility](http://alcf.anl.gov)
* [CALMIP](https://www.calmip.univ-toulouse.fr)
+
+------------------------------
+
+
+
+[TREX: Targeting Real Chemical Accuracy at the Exascale](https://trex-coe.eu) project has received funding from the European Union’s Horizon 2020 - Research and Innovation program - under grant agreement no. 952165. The content of this document does not represent the opinion of the European Union, and the European Union is not responsible for any use that might be made of such content.
+
+[PTEROSOR](https://lcpq.github.io/PTEROSOR) project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 863481).
diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio
index 07ad2236..e7c44b37 100755
--- a/bin/qp_convert_output_to_ezfio
+++ b/bin/qp_convert_output_to_ezfio
@@ -146,6 +146,17 @@ def write_ezfio(res, filename):
ezfio.set_ao_basis_ao_nucl(at)
ezfio.set_ao_basis_ao_prim_num(num_prim)
ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
+ try:
+ normf = res.normf
+ if normf == 0:
+ ezfio.set_ao_basis_ao_normalized(True)
+ elif normf == 1:
+ ezfio.set_ao_basis_ao_normalized(False)
+ else:
+ print("BUG in NORMF")
+ sys.exit(0)
+ except AttributeError:
+ ezfio.set_ao_basis_ao_normalized(True)
# ~#~#~#~#~#~#~ #
# P a r s i n g #
diff --git a/etc/qp.rc b/etc/qp.rc
index c56661c7..064ca3f7 100644
--- a/etc/qp.rc
+++ b/etc/qp.rc
@@ -80,8 +80,6 @@ function qp()
if [[ -d $NAME ]] ; then
[[ -d $EZFIO_FILE ]] && ezfio unset_file
ezfio set_file $NAME
- else
- qp_create_ezfio -h | more
fi
unset _ARGS
;;
diff --git a/ocaml/Command_line.ml b/ocaml/Command_line.ml
index 1dd57892..602315c6 100644
--- a/ocaml/Command_line.ml
+++ b/ocaml/Command_line.ml
@@ -1,3 +1,5 @@
+exception Error of string
+
type short_opt = char
type long_opt = string
type optional = Mandatory | Optional
@@ -181,15 +183,16 @@ let set_specs specs_in =
Getopt.parse_cmdline cmd_specs (fun x -> anon_args := !anon_args @ [x]);
if show_help () then
- (help () ; exit 0);
+ help ()
+ else
+ (* Check that all mandatory arguments are set *)
+ List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
+ |> List.iter (fun x ->
+ match get x.long with
+ | Some _ -> ()
+ | None -> raise (Error ("--"^x.long^" option is missing."))
+ )
- (* Check that all mandatory arguments are set *)
- List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs
- |> List.iter (fun x ->
- match get x.long with
- | Some _ -> ()
- | None -> failwith ("Error: --"^x.long^" option is missing.")
- )
;;
diff --git a/ocaml/Command_line.mli b/ocaml/Command_line.mli
index 9f6e7022..5ad4ee08 100644
--- a/ocaml/Command_line.mli
+++ b/ocaml/Command_line.mli
@@ -59,6 +59,8 @@ let () =
*)
+exception Error of string
+
type short_opt = char
type long_opt = string
diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml
index 9b01ac3a..603244c8 100644
--- a/ocaml/Molecule.ml
+++ b/ocaml/Molecule.ml
@@ -101,7 +101,7 @@ let to_string_general ~f m =
|> String.concat "\n"
let to_string =
- to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x)
+ to_string_general ~f:(fun x -> Atom.to_string ~units:Units.Angstrom x)
let to_xyz =
to_string_general ~f:Atom.to_xyz
@@ -113,7 +113,7 @@ let of_xyz_string
s =
let l = String_ext.split s ~on:'\n'
|> List.filter (fun x -> x <> "")
- |> list_map (fun x -> Atom.of_string units x)
+ |> list_map (fun x -> Atom.of_string ~units x)
in
let ne = ( get_charge {
nuclei=l ;
diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml
index be6c305b..4583b118 100644
--- a/ocaml/qp_create_ezfio.ml
+++ b/ocaml/qp_create_ezfio.ml
@@ -677,6 +677,7 @@ let run ?o b au c d m p cart xyz_file =
let () =
+ try (
let open Command_line in
begin
@@ -734,7 +735,7 @@ If a file with the same name as the basis set exists, this file will be read. O
let basis =
match Command_line.get "basis" with
- | None -> assert false
+ | None -> ""
| Some x -> x
in
@@ -773,10 +774,14 @@ If a file with the same name as the basis set exists, this file will be read. O
let xyz_filename =
match Command_line.anon_args () with
- | [x] -> x
- | _ -> (Command_line.help () ; failwith "input file is missing")
+ | [] -> failwith "input file is missing"
+ | x::_ -> x
in
run ?o:output basis au charge dummy multiplicity pseudo cart xyz_filename
+ )
+ with
+ | Failure txt -> Printf.eprintf "Fatal error: %s\n%!" txt
+ | Command_line.Error txt -> Printf.eprintf "Command line error: %s\n%!" txt
diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml
index d096b15b..dfbab167 100644
--- a/ocaml/qp_run.ml
+++ b/ocaml/qp_run.ml
@@ -110,7 +110,7 @@ let run slave ?prefix exe ezfio_file =
let task_thread =
let thread =
Thread.create ( fun () ->
- TaskServer.run port_number )
+ TaskServer.run ~port:port_number )
in
thread ();
in
diff --git a/src/cisd/cisd.irp.f b/src/cisd/cisd.irp.f
index 5f167686..3e1e8d97 100644
--- a/src/cisd/cisd.irp.f
+++ b/src/cisd/cisd.irp.f
@@ -47,6 +47,37 @@ program cisd
PROVIDE N_states
read_wf = .False.
SOFT_TOUCH read_wf
+
+ integer :: i,k
+
+ if(pseudo_sym)then
+ call H_apply_cisd_sym
+ else
+ call H_apply_cisd
+ endif
+ double precision :: r1, r2
+ double precision, allocatable :: U_csf(:,:)
+
+ allocate(U_csf(N_csf,N_states))
+ U_csf = 0.d0
+ U_csf(1,1) = 1.d0
+ do k=2,N_states
+ do i=1,N_csf
+ call random_number(r1)
+ call random_number(r2)
+ r1 = dsqrt(-2.d0*dlog(r1))
+ r2 = dacos(-1.d0)*2.d0*r2
+ U_csf(i,k) = r1*dcos(r2)
+ enddo
+ U_csf(k,k) = U_csf(k,k) +100.d0
+ enddo
+ do k=1,N_states
+ call normalize(U_csf(1,k),N_csf)
+ enddo
+ call convertWFfromCSFtoDET(N_states,U_csf(1,1),psi_coef(1,1))
+ deallocate(U_csf)
+ SOFT_TOUCH psi_coef
+
call run
end
@@ -56,13 +87,7 @@ subroutine run
double precision :: cisdq(N_states), delta_e
double precision,external :: diag_h_mat_elem
- if(pseudo_sym)then
- call H_apply_cisd_sym
- else
- call H_apply_cisd
- endif
psi_coef = ci_eigenvectors
- SOFT_TOUCH psi_coef
call save_wavefunction_truncated(save_threshold)
call ezfio_set_cisd_energy(CI_energy)
diff --git a/src/csf/cfgCI_utils.c b/src/csf/cfgCI_utils.c
index 23b984a0..746de04e 100644
--- a/src/csf/cfgCI_utils.c
+++ b/src/csf/cfgCI_utils.c
@@ -68,10 +68,16 @@ void getBFIndexList(int NSOMO, int *BF1, int *IdxListBF1){
break;
}
}
- BFcopy[Iidx] = -1;
- BFcopy[Jidx] = -1;
- IdxListBF1[Jidx] = Iidx;
- IdxListBF1[Iidx] = Jidx;
+ if(countN1 <= 0){
+ BFcopy[Iidx] = -1;
+ IdxListBF1[Iidx] = Iidx;
+ }
+ else{
+ BFcopy[Iidx] = -1;
+ BFcopy[Jidx] = -1;
+ IdxListBF1[Jidx] = Iidx;
+ IdxListBF1[Iidx] = Jidx;
+ }
}
}
@@ -328,10 +334,21 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
Get Overlap
************************************/
// Fill matrix
+
+ int rowsbftodetI, colsbftodetI;
+
+ /***********************************
+ Get BFtoDeterminant Matrix
+ ************************************/
+
+ printf("In convertcsftodet\n");
+ convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
+
int rowsI = 0;
int colsI = 0;
- getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+ //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+ getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
/***********************************
@@ -342,14 +359,6 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
gramSchmidt(overlapMatrixI, rowsI, colsI, orthoMatrixI);
- /***********************************
- Get BFtoDeterminant Matrix
- ************************************/
-
- int rowsbftodetI, colsbftodetI;
-
- convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
-
/***********************************
Get Final CSF to Det Matrix
************************************/
@@ -1297,12 +1306,16 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
double fac = 1.0;
for(int i = 0; i < NSOMO; i++)
donepq[i] = 0.0;
+ for(int i=0;i 0.0 || donepq[idxq] > 0.0) continue;
+ if(donepq[idxp] > 0.0 || donepq[idxq] > 0.0 || idxp == idxq) continue;
fac *= 2.0;
donepq[idxp] = 1.0;
donepq[idxq] = 1.0;
@@ -1319,15 +1332,19 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
}
shft /= 2;
}
-
+
// Now get the addresses
int inpdet[NSOMO];
int phase_cfg_to_qp=1;
int addr = -1;
for(int i = 0; i < npairs; i++){
- for(int j = 0; j < NSOMO; j++)
+ for(int j = 0; j < NSOMO; j++) {
inpdet[j] = detslist[i*NSOMO + j];
+ printf(" %d ",inpdet[j]);
+ }
+ printf("\n");
findAddofDetDriver(dettree, NSOMO, inpdet, &addr);
+ printf("(%d) - addr = %d\n",i,addr);
// Calculate the phase for cfg to QP2 conversion
//get_phase_cfg_to_qp_inpList(inpdet, NSOMO, &phase_cfg_to_qp);
//rowvec[addr] = 1.0 * phaselist[i]*phase_cfg_to_qp/sqrt(fac);
@@ -1416,7 +1433,6 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
getIthBFDriver(&bftree, NSOMO, addI, BF1);
getBFIndexList(NSOMO, BF1, IdxListBF1);
-
// Get ith row
getbftodetfunction(&dettree, NSOMO, MS, IdxListBF1, rowvec);
@@ -1425,6 +1441,11 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
for(int k=0;k NSOMOMax || icpl < 0 || izeros > zeromax ) return;
// If we find a valid BF assign its address
- if(isomo == NSOMOMax){
+ if(isomo == NSOMOMax && icpl == MSmax){
(*inode)->addr = bftree->rootNode->addr;
bftree->rootNode->addr += 1;
return;
}
// Call 0 branch
- if(((*inode)->C0) == NULL && izeros+1 <= zeromax){
- ((*inode)->C0) = malloc(sizeof(Node));
- (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
- buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
+ if(izeros+1 <= zeromax){
+ if(((*inode)->C0) == NULL){
+ ((*inode)->C0) = malloc(sizeof(Node));
+ (*(*inode)->C0) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 0, .iSOMO = isomo };
+ buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
+ }
+ else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
}
- else buildTree(bftree, &(*inode)->C0, isomo+1, izeros+1, icpl+1, NSOMOMax, MSmax);
// Call 1 branch
- if(((*inode)->C1) == NULL && icpl-1 >= 0){
- ((*inode)->C1) = malloc(sizeof(Node));
- (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
- buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
+ if(icpl-1 >=0){
+ if(((*inode)->C1) == NULL){
+ ((*inode)->C1) = malloc(sizeof(Node));
+ (*(*inode)->C1) = (Node){ .C0 = NULL, .C1 = NULL, .PREV = *inode, .addr = -1, .cpl = 1, .iSOMO = isomo };
+ buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
+ }
+ else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
}
- else buildTree(bftree, &(*inode)->C1, isomo+1, izeros+0, icpl-1, NSOMOMax, MSmax);
return;
}
diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f
index 0c3c6f92..7aaaa842 100644
--- a/src/davidson/diagonalization_hcsf_dressed.irp.f
+++ b/src/davidson/diagonalization_hcsf_dressed.irp.f
@@ -124,7 +124,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
stop -1
endif
- itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1
+ itermax = max(2,min(davidson_sze_max, sze_csf/N_st_diag))+1
itertot = 0
if (state_following) then
@@ -263,29 +263,20 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
! ===================
converged = .False.
-
+ call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),U_csf(1,1))
do k=N_st+1,N_st_diag
- do i=1,sze
+ do i=1,sze_csf
call random_number(r1)
call random_number(r2)
r1 = dsqrt(-2.d0*dlog(r1))
r2 = dtwo_pi*r2
- u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st)
+ U_csf(i,k) = r1*dcos(r2) * u_csf(i,k-N_st)
enddo
- u_in(k,k) = u_in(k,k) + 10.d0
+ U_csf(k,k) = u_csf(k,k) + 10.d0
enddo
do k=1,N_st_diag
- call normalize(u_in(1,k),sze)
+ call normalize(U_csf(1,k),sze_csf)
enddo
-
- do k=1,N_st_diag
- do i=1,sze
- U(i,k) = u_in(i,k)
- enddo
- enddo
-
- ! Make random verctors eigenstates of S2
- call convertWFfromDETtoCSF(N_st_diag,U(1,1),U_csf(1,1))
call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1))
do while (.not.converged)
diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f
index 46ad8f78..6930cc07 100644
--- a/src/davidson/diagonalize_ci.irp.f
+++ b/src/davidson/diagonalize_ci.irp.f
@@ -1,3 +1,13 @@
+BEGIN_PROVIDER [ character*(3), sigma_vector_algorithm ]
+ implicit none
+ BEGIN_DOC
+ ! If 'det', use in Davidson
+ !
+ ! If 'cfg', use in Davidson
+ END_DOC
+ sigma_vector_algorithm = 'det'
+END_PROVIDER
+
BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
implicit none
BEGIN_DOC
@@ -61,9 +71,18 @@ END_PROVIDER
if (diag_algorithm == "Davidson") then
if (do_csf) then
- call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
- size(CI_eigenvectors,1),CI_electronic_energy, &
- N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+ if (sigma_vector_algorithm == 'det') then
+ call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
+ size(CI_eigenvectors,1),CI_electronic_energy, &
+ N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+! else if (sigma_vector_algorithm == 'cfg') then
+! call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
+! size(CI_eigenvectors,1),CI_electronic_energy, &
+! N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
+! else
+! print *, irp_here
+! stop 'bug'
+ endif
else
call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_s2, &
size(CI_eigenvectors,1),CI_electronic_energy, &
diff --git a/src/determinants/determinants.irp.f b/src/determinants/determinants.irp.f
index 4b317025..e1c14bfe 100644
--- a/src/determinants/determinants.irp.f
+++ b/src/determinants/determinants.irp.f
@@ -77,7 +77,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
END_DOC
PROVIDE ezfio_filename
logical :: exists
- psi_det_size = 1
+ psi_det_size = N_states
PROVIDE mpi_master
if (read_wf) then
if (mpi_master) then
@@ -85,7 +85,7 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
if (exists) then
call ezfio_get_determinants_n_det(psi_det_size)
else
- psi_det_size = 1
+ psi_det_size = N_states
endif
call write_int(6,psi_det_size,'Dimension of the psi arrays')
endif
diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f
index 79897af7..56d8cf28 100644
--- a/src/mo_two_e_ints/mo_bi_integrals.irp.f
+++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f
@@ -130,7 +130,6 @@ subroutine four_idx_dgemm
real(integral_kind), allocatable :: buffer_value(:)
size_buffer = min(ao_num*ao_num*ao_num,16000000)
- print *, 'Storing'
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals)
allocate ( buffer_i(size_buffer), buffer_value(size_buffer) )
@@ -164,7 +163,6 @@ subroutine four_idx_dgemm
deallocate (a1)
- print *, 'Unique'
call map_unique(mo_integrals_map)
integer*8 :: get_mo_map_size, mo_map_size
diff --git a/src/tools/truncate_wf.irp.f b/src/tools/truncate_wf.irp.f
index 6c66c8ec..64c15bf7 100644
--- a/src/tools/truncate_wf.irp.f
+++ b/src/tools/truncate_wf.irp.f
@@ -54,11 +54,23 @@ subroutine routine_s2
double precision, allocatable :: psi_coef_tmp(:,:)
integer :: i,j,k
double precision :: accu(N_states)
+ integer :: weights(0:16), ix
+ double precision :: x
- print *, 'Weights of the CFG'
+ weights(:) = 0
do i=1,N_det
- print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
+ x = -dlog(1.d-32+sum(weight_configuration(det_to_configuration(i),:)))/dlog(10.d0)
+ ix = min(int(x), 16)
+ weights(ix) += 1
enddo
+
+ print *, 'Histogram of the weights of the CFG'
+ do i=0,15
+ print *, ' 10^{-', i, '} ', weights(i)
+ end do
+ print *, '< 10^{-', 15, '} ', weights(16)
+
+
print*, 'Min weight of the configuration?'
read(5,*) wmin