From ff86d51c5ffb195a6ebefa9d82599598233f7f53 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 25 Jun 2014 14:58:58 +0200 Subject: [PATCH] Hartree-Fock works well --- scripts/generate_h_apply.py | 24 +++-- setup_environment.sh | 4 +- src/AOs/ASSUMPTIONS.rst | 2 +- src/AOs/README.rst | 16 ++- src/BiInts/README.rst | 56 ++++++----- src/Bitmask/README.rst | 10 +- src/Bitmask/bitmasks.irp.f | 12 +++ src/CIS/ASSUMPTIONS.rst | 1 + src/CIS/H_apply.irp.f | 9 ++ src/CIS/Makefile | 8 ++ src/CIS/README.rst | 42 ++++++++ src/CIS/cis.irp.f | 17 ++++ src/CIS/super_ci.irp.f | 61 +++++++++++ src/DensityMatrix/README.rst | 6 -- src/Dets/README.rst | 18 ++++ src/Dets/density_matrix.irp.f | 15 ++- src/Dets/determinants.ezfio_config | 4 + src/Dets/determinants.irp.f | 2 +- src/Dets/save_for_qmcchem.irp.f | 37 +++++++ src/Generators_full/README.rst | 42 ++++++++ src/Generators_full/generators.irp.f | 10 +- src/Hartree_Fock/Fock_matrix.irp.f | 33 ++++++ src/Hartree_Fock/README.rst | 28 +----- src/Hartree_Fock/SCF.irp.f | 117 +++------------------ src/Hartree_Fock/damping_SCF.irp.f | 123 +++++++++++++++++++++++ src/Hartree_Fock/mo_SCF_iterations.irp.f | 23 ----- src/MOGuess/NEEDED_MODULES | 3 +- src/MOGuess/README.rst | 11 +- src/MOs/README.rst | 5 +- src/MOs/mos.ezfio_config | 1 - src/Makefile.config.example | 2 +- src/Output/output.irp.f | 2 +- src/Perturbation/README.rst | 20 ++-- src/SC2/README.rst | 58 ----------- src/Selectors_full/README.rst | 93 +++++++++++++++++ src/SingleRefMethod/README.rst | 3 + src/SingleRefMethod/generators.irp.f | 1 + src/Utils/LinearAlgebra.irp.f | 122 +++++++++++++++++++++- src/Utils/README.rst | 6 +- 39 files changed, 751 insertions(+), 296 deletions(-) create mode 100644 src/CIS/ASSUMPTIONS.rst create mode 100644 src/CIS/H_apply.irp.f create mode 100644 src/CIS/Makefile create mode 100644 src/CIS/README.rst create mode 100644 src/CIS/cis.irp.f create mode 100644 src/CIS/super_ci.irp.f create mode 100644 src/Dets/save_for_qmcchem.irp.f create mode 100644 src/Hartree_Fock/damping_SCF.irp.f delete mode 100644 src/Hartree_Fock/mo_SCF_iterations.irp.f diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index a57f869f..083a91ba 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -216,8 +216,9 @@ class H_apply(object): self.data["size_max"] = str(1024*128) self.data["copy_buffer"] = """ call copy_h_apply_buffer_to_wf - selection_criterion_min = selection_criterion_min*0.1d0 + selection_criterion_min = min(selection_criterion_min, maxval(select_max))*0.1d0 selection_criterion = selection_criterion_min + call write_double(output_Dets,selection_criterion,'Selection criterion') """ self.data["keys_work"] = """ e_2_pert_buffer = 0.d0 @@ -230,16 +231,17 @@ class H_apply(object): self.data["omp_parallel"] += """& !$OMP REDUCTION (max:select_max_out)""" self.data["skip"] = """ - if ((i_generator < size(select_max)).and. & - (select_max(i_generator) < selection_criterion_min*selection_criterion_factor)) then - !$ call omp_set_lock(lck) - do k=1,N_st - norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k) - delta_pt2(k) = 0.d0 - pt2_old(k) = 0.d0 - enddo - !$ call omp_unset_lock(lck) - cycle + if (i_generator < size_select_max) then + if (select_max(i_generator) < selection_criterion_min*selection_criterion_factor) then + !$ call omp_set_lock(lck) + do k=1,N_st + norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k) + delta_pt2(k) = 0.d0 + pt2_old(k) = 0.d0 + enddo + !$ call omp_unset_lock(lck) + cycle + endif endif select_max(i_generator) = 0.d0 """ diff --git a/setup_environment.sh b/setup_environment.sh index 253a1043..76d82885 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -23,9 +23,9 @@ export QPACKAGE_ROOT=${QPACKAGE_ROOT} export PYTHONPATH+=:\${QPACKAGE_ROOT}/scripts export PATH+=:\${QPACKAGE_ROOT}/scripts export PATH+=:\${QPACKAGE_ROOT}/bin -export QPACKAGE_CACHE_URL="http://qmcchem.ups-tlse.fr/files/scemama/quantum_package/cache +export QPACKAGE_CACHE_URL="http://qmcchem.ups-tlse.fr/files/scemama/quantum_package/cache" export PATH+=:${QPACKAGE_ROOT}/irpf90/bin/ -" +source ${QPACKAGE_ROOT}/irpf90/bin/irpman EOF source quantum_package.rc diff --git a/src/AOs/ASSUMPTIONS.rst b/src/AOs/ASSUMPTIONS.rst index cb89b682..0c8e0ccc 100644 --- a/src/AOs/ASSUMPTIONS.rst +++ b/src/AOs/ASSUMPTIONS.rst @@ -4,5 +4,5 @@ \int \left(\chi_i({\bf r}) \right)^2 dr = 1 -* The AO coefficients in the EZFIO files are not normalized +* The AO coefficients in the EZFIO files are not necessarily normalized and are normalized after reading * The AO coefficients and exponents are ordered in increasing order of exponents diff --git a/src/AOs/README.rst b/src/AOs/README.rst index 378bf8ae..e01108cd 100644 --- a/src/AOs/README.rst +++ b/src/AOs/README.rst @@ -54,17 +54,17 @@ Documentation Coefficients, exponents and powers of x,y and z ao_coef(i,j) = coefficient of the jth primitive on the ith ao -`ao_coef_transp `_ +`ao_coef_transp `_ Transposed ao_coef and ao_expo `ao_expo `_ Coefficients, exponents and powers of x,y and z ao_coef(i,j) = coefficient of the jth primitive on the ith ao -`ao_expo_transp `_ +`ao_expo_transp `_ Transposed ao_coef and ao_expo -`ao_nucl `_ +`ao_nucl `_ Index of the nuclei on which the ao is centered `ao_num `_ @@ -73,21 +73,17 @@ Documentation `ao_num_align `_ Number of atomic orbitals -`ao_overlap `_ - matrix of the overlap for tha AOs - S(i,j) = overlap between the ith and the jth atomic basis function - `ao_power `_ Coefficients, exponents and powers of x,y and z ao_coef(i,j) = coefficient of the jth primitive on the ith ao -`ao_prim_num `_ +`ao_prim_num `_ Number of primitives per atomic orbital -`ao_prim_num_max `_ +`ao_prim_num_max `_ Undocumented -`ao_prim_num_max_align `_ +`ao_prim_num_max_align `_ Undocumented diff --git a/src/BiInts/README.rst b/src/BiInts/README.rst index 48fa022b..e71ab5ba 100644 --- a/src/BiInts/README.rst +++ b/src/BiInts/README.rst @@ -36,58 +36,58 @@ Documentation integral of the AO basis or (ij|kl) i(r1) j(r1) 1/r12 k(r2) l(r2) -`ao_bielec_integral_schwartz `_ +`ao_bielec_integral_schwartz `_ Needed to compuet Schwartz inequalities -`ao_bielec_integrals_in_map `_ +`ao_bielec_integrals_in_map `_ Map of Atomic integrals i(r1) j(r2) 1/r12 k(r1) l(r2) -`compute_ao_bielec_integrals `_ +`compute_ao_bielec_integrals `_ Compute AO 1/r12 integrals for all i and fixed j,k,l -`eri `_ +`eri `_ ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) -`general_primitive_integral `_ +`general_primitive_integral `_ Computes the integral where p,q,r,s are Gaussian primitives -`give_polynom_mult_center_x `_ +`give_polynom_mult_center_x `_ subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw : I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) -`i_x1_new `_ +`i_x1_new `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult `_ +`i_x1_pol_mult `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a1 `_ +`i_x1_pol_mult_a1 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a2 `_ +`i_x1_pol_mult_a2 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_recurs `_ +`i_x1_pol_mult_recurs `_ recursive function involved in the bielectronic integral -`i_x2_new `_ +`i_x2_new `_ recursive function involved in the bielectronic integral -`i_x2_pol_mult `_ +`i_x2_pol_mult `_ recursive function involved in the bielectronic integral -`integrale_new `_ +`integrale_new `_ calculate the integral of the polynom :: I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q) between ( 0 ; 1) -`n_pt_sup `_ +`n_pt_sup `_ Returns the upper boundary of the degree of the polynom involved in the bielctronic integral : Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z) @@ -112,7 +112,7 @@ Documentation `clear_ao_map `_ Frees the memory of the AO map -`clear_mo_map `_ +`clear_mo_map `_ Frees the memory of the MO map `get_ao_bielec_integral `_ @@ -136,7 +136,11 @@ Documentation Returns multiple integrals in the MO basis, all i for j,k,l fixed. -`get_mo_map_size `_ +`get_mo_bielec_integrals_existing_ik `_ + Returns multiple integrals in the MO basis, all + i for j,k,l fixed. + +`get_mo_map_size `_ Return the number of elements in the MO map `insert_into_ao_integrals_map `_ @@ -154,13 +158,13 @@ Documentation `add_integrals_to_map `_ Adds integrals to tha MO map according to some bitmask -`mo_bielec_integral_jj `_ +`mo_bielec_integral_jj `_ Transform Bi-electronic integrals and -`mo_bielec_integral_jj_anti `_ +`mo_bielec_integral_jj_anti `_ Transform Bi-electronic integrals and -`mo_bielec_integral_jj_exchange `_ +`mo_bielec_integral_jj_exchange `_ Transform Bi-electronic integrals and `mo_bielec_integrals_in_map `_ @@ -169,22 +173,22 @@ Documentation `mo_bielec_integrals_index `_ Computes an unique index for i,j,k,l integrals -`ao_integrals_threshold `_ +`ao_integrals_threshold `_ If < ao_integrals_threshold, = 0 -`do_direct_integrals `_ +`do_direct_integrals `_ If True, compute integrals on the fly -`mo_integrals_threshold `_ +`mo_integrals_threshold `_ If < mo_integrals_threshold, = 0 -`read_ao_integrals `_ +`read_ao_integrals `_ If true, read AO integrals in EZFIO -`read_mo_integrals `_ +`read_mo_integrals `_ If true, read MO integrals in EZFIO -`write_ao_integrals `_ +`write_ao_integrals `_ If true, write AO integrals in EZFIO `write_mo_integrals `_ diff --git a/src/Bitmask/README.rst b/src/Bitmask/README.rst index 057bd7eb..01b198ab 100644 --- a/src/Bitmask/README.rst +++ b/src/Bitmask/README.rst @@ -57,7 +57,7 @@ Documentation `full_ijkl_bitmask `_ Bitmask to include all possible MOs -`generators_bitmask `_ +`generators_bitmask `_ Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator). 3rd index is : * 1 : hole for single exc @@ -70,10 +70,10 @@ Documentation `hf_bitmask `_ Hartree Fock bit mask -`i_bitmask_gen `_ +`i_bitmask_gen `_ Current bitmask for the generators -`i_bitmask_ref `_ +`i_bitmask_ref `_ Current bitmask for the reference `n_generators_bitmask `_ @@ -82,13 +82,13 @@ Documentation `n_int `_ Number of 64-bit integers needed to represent determinants as binary strings -`n_reference_bitmask `_ +`n_reference_bitmask `_ Number of bitmasks for reference `ref_bitmask `_ Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask -`reference_bitmask `_ +`reference_bitmask `_ Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference) `bitstring_to_hexa `_ diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index bec60e53..2f841a7d 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -28,6 +28,18 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask, (N_int,4) ] enddo END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), cis_ijkl_bitmask, (N_int,4) ] + implicit none + BEGIN_DOC + ! Bitmask to include all possible single excitations from Hartree-Fock + END_DOC + + integer :: i,j,n + cis_ijkl_bitmask = full_ijkl_bitmask + cis_ijkl_bitmask(:,1) = HF_bitmask(:,1) +END_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask, (N_int,2)] implicit none diff --git a/src/CIS/ASSUMPTIONS.rst b/src/CIS/ASSUMPTIONS.rst new file mode 100644 index 00000000..bbe7afec --- /dev/null +++ b/src/CIS/ASSUMPTIONS.rst @@ -0,0 +1 @@ +* The molecular orbitals are assumed orthonormal diff --git a/src/CIS/H_apply.irp.f b/src/CIS/H_apply.irp.f new file mode 100644 index 00000000..cf68267e --- /dev/null +++ b/src/CIS/H_apply.irp.f @@ -0,0 +1,9 @@ +! Generates subroutine H_apply_cisd +! ---------------------------------- + +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import H_apply +H = H_apply("cis",do_double_exc=False) +print H +END_SHELL + diff --git a/src/CIS/Makefile b/src/CIS/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/CIS/Makefile @@ -0,0 +1,8 @@ +default: all + +# Define here all new external source files and objects.Don't forget to prefix the +# object files with IRPF90_temp/ +SRC= +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/CIS/README.rst b/src/CIS/README.rst new file mode 100644 index 00000000..47e768c9 --- /dev/null +++ b/src/CIS/README.rst @@ -0,0 +1,42 @@ +========== +CIS Module +========== + +Assumptions +=========== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* The molecular orbitals are assumed orthonormal + + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `SingleRefMethod `_ +* `Utils `_ +* `Selectors_full `_ + diff --git a/src/CIS/cis.irp.f b/src/CIS/cis.irp.f new file mode 100644 index 00000000..4859c0d1 --- /dev/null +++ b/src/CIS/cis.irp.f @@ -0,0 +1,17 @@ +program cis + implicit none + integer :: i + + print *, 'HF = ', HF_energy + print *, 'N_states = ', N_states + call H_apply_cis + print *, 'N_det = ', N_det + do i = 1,N_states + print *, 'energy = ',CI_energy(i) + print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy + enddo + psi_coef = ci_eigenvectors + SOFT_TOUCH psi_coef + call save_wavefunction + +end diff --git a/src/CIS/super_ci.irp.f b/src/CIS/super_ci.irp.f new file mode 100644 index 00000000..a0b33c14 --- /dev/null +++ b/src/CIS/super_ci.irp.f @@ -0,0 +1,61 @@ +program cis + implicit none + integer :: i + + call super_CI + +end + +subroutine super_CI + implicit none + double precision :: E, delta_E, delta_D, E_min + integer :: k + character :: save_char + + call write_time(output_hartree_fock) + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'), & + '====','================','================','================' + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'), & + ' N ', 'Energy ', 'Energy diff ', 'Save ' + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'), & + '====','================','================','================' + + E = HF_energy + 1.d0 + delta_D = 0.d0 + E_min = HF_energy + FREE psi_det psi_coef + call clear_mo_map + N_det = 1 + SOFT_TOUCH N_det + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + do k=1,n_it_scf_max + delta_E = HF_energy - E + E = HF_energy + if (E < E_min) then + call save_mos + save_char = 'X' + else + save_char = ' ' + endif + write(output_hartree_fock,'(I4,X,F16.10, X, F16.10, X, A8 )'),& + k, E, delta_E, save_char + if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then + exit + endif + call H_apply_cis + call diagonalize_CI + call set_natural_mos + FREE psi_det psi_coef + call clear_mo_map + N_det = 1 + SOFT_TOUCH N_det + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + enddo + + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16 )'), & + '====','================','================','================' + call write_time(output_hartree_fock) +end + diff --git a/src/DensityMatrix/README.rst b/src/DensityMatrix/README.rst index f232ad4e..a8c58436 100644 --- a/src/DensityMatrix/README.rst +++ b/src/DensityMatrix/README.rst @@ -17,12 +17,6 @@ Documentation `iunit_two_body_dm_bb `_ Temporary files for 2-body dm calculation -`one_body_dm_a `_ - Alpha and beta one-body density matrix - -`one_body_dm_b `_ - Alpha and beta one-body density matrix - `two_body_dm_diag_aa `_ diagonal part of the two body density matrix diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 4a4e8ff2..4f9cda62 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -152,6 +152,21 @@ Documentation `davidson_threshold `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] +`one_body_dm_mo `_ + One-body density matrix + +`one_body_dm_mo_alpha `_ + Alpha and beta one-body density matrix for each state + +`one_body_dm_mo_beta `_ + Alpha and beta one-body density matrix for each state + +`save_natural_mos `_ + Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + +`state_average_weight `_ + Weights in the state-average calculation of the density matrix + `det_search_key `_ Return an integer*8 corresponding to a determinant index for searching @@ -307,6 +322,9 @@ Documentation `s_z2_sz `_ Undocumented +`save_natorb `_ + Undocumented + `a_operator `_ Needed for diag_H_mat_elem diff --git a/src/Dets/density_matrix.irp.f b/src/Dets/density_matrix.irp.f index e5a907fb..8fe608ab 100644 --- a/src/Dets/density_matrix.irp.f +++ b/src/Dets/density_matrix.irp.f @@ -78,10 +78,10 @@ BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta END_PROVIDER -subroutine save_natural_mos +subroutine set_natural_mos implicit none BEGIN_DOC - ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + ! Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis END_DOC character*(64) :: label double precision, allocatable :: tmp(:,:) @@ -92,9 +92,18 @@ subroutine save_natural_mos label = "Natural" call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label) deallocate(tmp) - call save_mos end +subroutine save_natural_mos + implicit none + BEGIN_DOC + ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + END_DOC + call set_natural_mos + call save_mos + +end + BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] implicit none diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index 5a01925a..b43a2dcb 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -10,3 +10,7 @@ determinants n_det_max_jacobi integer threshold_generators double precision threshold_selectors double precision + det_num integer + det_occ integer (electrons_elec_alpha_num,determinants_det_num,2) + det_coef double precision (determinants_det_num) + diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index d93507c0..e076267d 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -56,7 +56,7 @@ BEGIN_PROVIDER [ integer, N_det_max_jacobi ] if (exists) then call ezfio_get_determinants_n_det_max_jacobi(N_det_max_jacobi) else - N_det_max_jacobi = 2000 + N_det_max_jacobi = 10000 endif call write_int(output_dets,N_det_max_jacobi,'Lapack diagonalization up to') ASSERT (N_det_max_jacobi > 0) diff --git a/src/Dets/save_for_qmcchem.irp.f b/src/Dets/save_for_qmcchem.irp.f new file mode 100644 index 00000000..b6f60eb6 --- /dev/null +++ b/src/Dets/save_for_qmcchem.irp.f @@ -0,0 +1,37 @@ +subroutine save_dets_qmcchem + use bitmasks + implicit none + character :: c(mo_tot_num) + integer :: i,k + + integer, allocatable :: occ(:,:,:), occ_tmp(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: occ, occ_tmp + + call ezfio_set_determinants_det_num(N_det) + call ezfio_set_determinants_det_coef(psi_coef_sorted(1,1)) + + allocate (occ(elec_alpha_num,N_det,2)) + ! OMP PARALLEL DEFAULT(NONE) & + ! OMP PRIVATE(occ_tmp,i,k)& + ! OMP SHARED(N_det,psi_det_sorted,elec_alpha_num, & + ! OMP occ,elec_beta_num,N_int) + allocate (occ_tmp(N_int*bit_kind_size,2)) + occ_tmp = 0 + ! OMP DO + do i=1,N_det + call bitstring_to_list(psi_det_sorted(1,1,i), occ_tmp(1,1), elec_alpha_num, N_int ) + call bitstring_to_list(psi_det_sorted(1,2,i), occ_tmp(1,2), elec_beta_num, N_int ) + do k=1,elec_alpha_num + occ(k,i,1) = occ_tmp(k,1) + occ(k,i,2) = occ_tmp(k,2) + enddo + enddo + ! OMP END DO + deallocate(occ_tmp) + ! OMP END PARALLEL + call ezfio_set_determinants_det_occ(occ) + call write_int(output_dets,N_det,'Determinants saved for QMC') + deallocate(occ) +end + + diff --git a/src/Generators_full/README.rst b/src/Generators_full/README.rst index e3ba2312..3e7f8134 100644 --- a/src/Generators_full/README.rst +++ b/src/Generators_full/README.rst @@ -5,3 +5,45 @@ Generators_full Module All the determinants of the wave function are generators. In this way, the Full CI space is explored. +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`n_det_generators `_ + For Single reference wave functions, the number of generators is 1 : the + Hartree-Fock determinant + +`psi_generators `_ + For Single reference wave functions, the generator is the + Hartree-Fock determinant + +`select_max `_ + Memo to skip useless selectors + +`threshold_generators `_ + Percentage of the norm of the state-averaged wave function to + consider for the generators + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Utils `_ + diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index 9d8761d2..e0db1fc3 100644 --- a/src/Generators_full/generators.irp.f +++ b/src/Generators_full/generators.irp.f @@ -57,7 +57,15 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] END_PROVIDER - BEGIN_PROVIDER [ double precision, select_max, (3000) ] +BEGIN_PROVIDER [ integer, size_select_max] + implicit none + BEGIN_DOC + ! Size of the select_max array + END_DOC + size_select_max = 10000 +END_PROVIDER + +BEGIN_PROVIDER [ double precision precision, select_max, (size_select_max) ] implicit none BEGIN_DOC ! Memo to skip useless selectors diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 69aa7e26..3f39ee24 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -260,3 +260,36 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] endif END_PROVIDER +subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) + implicit none + integer, intent(in) :: LDFMO ! size(FMO,1) + integer, intent(in) :: LDFAO ! size(FAO,1) + double precision, intent(in) :: FMO(LDFMO,*) + double precision, intent(out) :: FAO(LDFAO,*) + + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + FMO, size(FMO,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + FAO, size(FAO,1)) + deallocate(T,M) +end + diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 341c6284..3d77b7de 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -71,6 +71,9 @@ Documentation K = Fb - Fa .br +`fock_mo_to_ao `_ + Undocumented + `hf_energy `_ Hartree-Fock energy @@ -83,25 +86,10 @@ Documentation `hf_density_matrix_ao_beta `_ Beta density matrix in the AO basis -`fock_mo_to_ao `_ +`run `_ Undocumented -`insert_new_scf_density_matrix `_ - Undocumented - -`it_scf `_ - Number of the current SCF iteration - -`scf_density_matrices `_ - Density matrices at every SCF iteration - -`scf_energies `_ - Density matrices at every SCF iteration - -`scf_interpolation_step `_ - Undocumented - -`scf_iterations `_ +`scf `_ Undocumented `diagonal_fock_matrix_mo `_ @@ -110,12 +98,6 @@ Documentation `eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`run `_ - Undocumented - -`scf `_ - Undocumented - `do_diis `_ If True, compute integrals on the fly diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f index 650aa365..55c42178 100644 --- a/src/Hartree_Fock/SCF.irp.f +++ b/src/Hartree_Fock/SCF.irp.f @@ -1,108 +1,23 @@ -BEGIN_PROVIDER [ integer, it_scf ] - implicit none - BEGIN_DOC - ! Number of the current SCF iteration - END_DOC - it_scf = 0 -END_PROVIDER - BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ] -&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ] - implicit none - BEGIN_DOC - ! Density matrices at every SCF iteration - END_DOC - SCF_density_matrices = 0.d0 - SCF_energies = 0.d0 -END_PROVIDER - -subroutine insert_new_SCF_density_matrix - implicit none - integer :: i,j - do j=1,ao_num - do i=1,ao_num - SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) - SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) - enddo - enddo - SCF_energies(it_scf) = HF_energy +program scf + call orthonormalize_mos + call run end -subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) +subroutine run + + use bitmasks implicit none - integer, intent(in) :: LDFMO ! size(FMO,1) - integer, intent(in) :: LDFAO ! size(FAO,1) - double precision, intent(in) :: FMO(LDFMO,*) - double precision, intent(out) :: FAO(LDFAO,*) + double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral + double precision :: E0 + integer :: i_it, i, j, k + + E0 = HF_energy - double precision, allocatable :: T(:,:), M(:,:) - ! F_ao = S C F_mo C^t S - allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - ao_overlap, size(ao_overlap,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & - M, size(M,1), & - FMO, size(FMO,1), & - 0.d0, & - T, size(T,1)) - call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & - T, size(T,1), & - mo_coef, size(mo_coef,1), & - 0.d0, & - M, size(M,1)) - call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & - M, size(M,1), & - ao_overlap, size(ao_overlap,1), & - 0.d0, & - FAO, size(FAO,1)) - deallocate(T,M) -end - -subroutine SCF_interpolation_step - implicit none - integer :: i,j - double precision :: c + thresh_SCF = 1.d-10 + call damping_SCF + mo_label = "Canonical" + TOUCH mo_label mo_coef + call save_mos - if (it_scf == 1) then - return - endif - call random_number(c) - c = c*0.5d0 - do j=1,ao_num - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf-1)+SCF_density_matrices(i,j,1,it_scf) * (1.d0 - c) - HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf-1)+SCF_density_matrices(i,j,2,it_scf) * (1.d0 - c) - enddo - enddo - TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta - - ! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1),& - ! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) - ! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1),& - ! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) - ! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta -end - -subroutine scf_iterations - implicit none - integer :: i,j - do i=1,n_it_scf_max - it_scf += 1 - SOFT_TOUCH it_scf - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef - call insert_new_SCF_density_matrix - print *, HF_energy - if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then - call SCF_interpolation_step - endif - if (it_scf>1 ) then - if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then - exit - endif - endif - enddo end diff --git a/src/Hartree_Fock/damping_SCF.irp.f b/src/Hartree_Fock/damping_SCF.irp.f new file mode 100644 index 00000000..0e2f8f59 --- /dev/null +++ b/src/Hartree_Fock/damping_SCF.irp.f @@ -0,0 +1,123 @@ +subroutine damping_SCF + implicit none + double precision :: E + double precision, allocatable :: D_alpha(:,:), D_beta(:,:) + double precision :: E_new + double precision, allocatable :: D_new_alpha(:,:), D_new_beta(:,:), F_new(:,:) + double precision, allocatable :: delta_alpha(:,:), delta_beta(:,:) + double precision :: lambda, E_half, a, b, delta_D, delta_E, E_min + + integer :: i,j,k + logical :: saving + character :: save_char + + allocate( & + D_alpha( ao_num_align, ao_num ), & + D_beta( ao_num_align, ao_num ), & + F_new( ao_num_align, ao_num ), & + D_new_alpha( ao_num_align, ao_num ), & + D_new_beta( ao_num_align, ao_num ), & + delta_alpha( ao_num_align, ao_num ), & + delta_beta( ao_num_align, ao_num )) + + do j=1,ao_num + do i=1,ao_num + D_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) + D_beta (i,j) = HF_density_matrix_ao_beta (i,j) + enddo + enddo + + + call write_time(output_hartree_fock) + + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '====' + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), ' N ', 'Energy ', 'Energy diff ', 'Density diff ', 'Save' + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '====' + + E = HF_energy + 1.d0 + E_min = HF_energy + delta_D = 0.d0 + do k=1,n_it_scf_max + + delta_E = HF_energy - E + E = HF_energy + + if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then + exit + endif + + saving = E < E_min + if (saving) then + call save_mos + save_char = 'X' + E_min = E + else + save_char = ' ' + endif + + write(output_hartree_fock,'(I4,X,F16.10, X, F16.10, X, F16.10, 3X, A )'), & + k, E, delta_E, delta_D, save_char + + D_alpha = HF_density_matrix_ao_alpha + D_beta = HF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + + D_new_alpha = HF_density_matrix_ao_alpha + D_new_beta = HF_density_matrix_ao_beta + F_new = Fock_matrix_ao + E_new = HF_energy + + delta_alpha = D_new_alpha - D_alpha + delta_beta = D_new_beta - D_beta + + lambda = 0.5d0 + E_half = 0.d0 + do while (E_half > E) + HF_density_matrix_ao_alpha = D_alpha + lambda * delta_alpha + HF_density_matrix_ao_beta = D_beta + lambda * delta_beta + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + E_half = HF_energy + if ((E_half > E).and.(E_new < E)) then + lambda = 1.d0 + exit + else if ((E_half > E).and.(lambda > 1.d-3)) then + lambda = 0.5d0 * lambda + E_new = E_half + else + exit + endif + enddo + + a = (E_new + E - 2.d0*E_half)*2.d0 + b = -E_new - 3.d0*E + 4.d0*E_half + lambda = -lambda*b/a + D_alpha = (1.d0-lambda) * D_alpha + lambda * D_new_alpha + D_beta = (1.d0-lambda) * D_beta + lambda * D_new_beta + delta_E = HF_energy - E + do j=1,ao_num + do i=1,ao_num + delta_D = delta_D + & + (D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j))*(D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j)) + & + (D_beta (i,j) - HF_density_matrix_ao_beta (i,j))*(D_beta (i,j) - HF_density_matrix_ao_beta (i,j)) + enddo + enddo + delta_D = dsqrt(delta_D/dble(ao_num)**2) + HF_density_matrix_ao_alpha = D_alpha + HF_density_matrix_ao_beta = D_beta + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + mo_coef = eigenvectors_fock_matrix_mo + TOUCH mo_coef + + + enddo + write(output_hartree_fock,'(A4,X,A16, X, A16, X, A16, X, A4 )'), '====','================','================','================', '====' + write(output_hartree_fock,*) + + call write_double(output_hartree_fock, HF_energy, 'Hartree-Fock energy') + call write_time(output_hartree_fock) + + deallocate(D_alpha,D_beta,F_new,D_new_alpha,D_new_beta,delta_alpha,delta_beta) +end diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f deleted file mode 100644 index 55c42178..00000000 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ /dev/null @@ -1,23 +0,0 @@ - -program scf - call orthonormalize_mos - call run -end - -subroutine run - - use bitmasks - implicit none - double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral - double precision :: E0 - integer :: i_it, i, j, k - - E0 = HF_energy - - thresh_SCF = 1.d-10 - call damping_SCF - mo_label = "Canonical" - TOUCH mo_label mo_coef - call save_mos - -end diff --git a/src/MOGuess/NEEDED_MODULES b/src/MOGuess/NEEDED_MODULES index 27c50e45..b4266ae1 100644 --- a/src/MOGuess/NEEDED_MODULES +++ b/src/MOGuess/NEEDED_MODULES @@ -1,2 +1,3 @@ -AOs Ezfio_files MOs Nuclei Output Utils MonoInts +AOs Electrons Ezfio_files MonoInts MOs Nuclei Output Utils + diff --git a/src/MOGuess/README.rst b/src/MOGuess/README.rst index 90785358..cb1702ab 100644 --- a/src/MOGuess/README.rst +++ b/src/MOGuess/README.rst @@ -9,12 +9,13 @@ Needed Modules .. NEEDED_MODULES file. * `AOs `_ +* `Electrons `_ * `Ezfio_files `_ +* `MonoInts `_ * `MOs `_ * `Nuclei `_ * `Output `_ * `Utils `_ -* `MonoInts `_ Documentation ============= @@ -22,19 +23,19 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`h_core_guess `_ +`h_core_guess `_ Undocumented -`ao_ortho_lowdin_coef `_ +`ao_ortho_lowdin_coef `_ matrix of the coefficients of the mos generated by the orthonormalization by the S^{-1/2} canonical transformation of the aos ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital -`ao_ortho_lowdin_overlap `_ +`ao_ortho_lowdin_overlap `_ overlap matrix of the ao_ortho_lowdin supposed to be the Identity -`ao_ortho_lowdin_nucl_elec_integral `_ +`ao_ortho_lowdin_nucl_elec_integral `_ Undocumented diff --git a/src/MOs/README.rst b/src/MOs/README.rst index f27425bc..9e3e7e92 100644 --- a/src/MOs/README.rst +++ b/src/MOs/README.rst @@ -36,12 +36,15 @@ Documentation .. NEEDED_MODULES file. `cholesky_mo `_ - Cholesky decomposition of MO Density matrix to + Cholesky decomposition of AO Density matrix to generate MOs `mo_density_matrix `_ Density matrix in MO basis +`mo_density_matrix_virtual `_ + Density matrix in MO basis (virtual MOs) + `mo_coef `_ Molecular orbital coefficients on AO basis set mo_coef(i,j) = coefficient of the ith ao on the jth mo diff --git a/src/MOs/mos.ezfio_config b/src/MOs/mos.ezfio_config index e34d1aab..b43f23b9 100644 --- a/src/MOs/mos.ezfio_config +++ b/src/MOs/mos.ezfio_config @@ -1,5 +1,4 @@ mo_basis - ao_num integer mo_tot_num integer mo_coef double precision (ao_basis_ao_num,mo_basis_mo_tot_num) mo_label character*(64) diff --git a/src/Makefile.config.example b/src/Makefile.config.example index d576e98c..ef7c48c0 100644 --- a/src/Makefile.config.example +++ b/src/Makefile.config.example @@ -1,6 +1,6 @@ OPENMP =1 PROFILE =0 -DEBUG = 1 +DEBUG = 0 IRPF90_FLAGS+= --align=32 FC = ifort -g diff --git a/src/Output/output.irp.f b/src/Output/output.irp.f index 9218bc0a..1a025c73 100644 --- a/src/Output/output.irp.f +++ b/src/Output/output.irp.f @@ -55,7 +55,7 @@ subroutine write_double(iunit,value,label) integer, intent(in) :: iunit double precision :: value character*(*) :: label - character*(64), parameter :: f = '(A50,G16.8)' + character*(64), parameter :: f = '(A50,G24.16)' character*(50) :: newlabel write(newlabel,'(A,A)') '* ',trim(label) write(iunit,f) newlabel, value diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index aad1d0c0..75668b37 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -82,7 +82,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`pt2_moller_plesset `_ +`pt2_moller_plesset `_ compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution .br for the various n_st states. @@ -92,7 +92,7 @@ Documentation e_2_pert(i) = ^2/(difference of orbital energies) .br -`pt2_epstein_nesbet `_ +`pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states. @@ -102,7 +102,7 @@ Documentation e_2_pert(i) = ^2/( E(i) - ) .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various N_st states. @@ -112,7 +112,7 @@ Documentation c_pert(i) = e_2_pert(i)/ .br -`pt2_epstein_nesbet_sc2_projected `_ +`pt2_epstein_nesbet_sc2_projected `_ compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states, @@ -135,23 +135,23 @@ Documentation .br H_pert_diag = c_pert -`repeat_all_e_corr `_ +`repeat_all_e_corr `_ Undocumented -`fill_h_apply_buffer_selection `_ +`fill_h_apply_buffer_selection `_ Fill the H_apply buffer with determiants for the selection -`remove_small_contributions `_ +`remove_small_contributions `_ Remove determinants with small contributions. N_states is assumed to be provided. -`selection_criterion `_ +`selection_criterion `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_factor `_ +`selection_criterion_factor `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_min `_ +`selection_criterion_min `_ Threshold to select determinants. Set by selection routines. diff --git a/src/SC2/README.rst b/src/SC2/README.rst index 30269d2d..32a5c02f 100644 --- a/src/SC2/README.rst +++ b/src/SC2/README.rst @@ -8,67 +8,9 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd_sc2 `_ - CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) - .br - dets_in : bitmasks corresponding to determinants - .br - u_in : guess coefficients on the various states. Overwritten - on exit - .br - dim_in : leftmost dimension of u_in - .br - sze : Number of determinants - .br - N_st : Number of eigenstates - .br - Initial guess vectors are not necessarily orthonormal - -`repeat_excitation `_ - Undocumented - `cisd `_ Undocumented -`ci_sc2_eigenvectors `_ - Eigenvectors/values of the CI matrix - -`ci_sc2_electronic_energy `_ - Eigenvectors/values of the CI matrix - -`ci_sc2_energy `_ - N_states lowest eigenvalues of the CI matrix - -`diagonalize_ci_sc2 `_ - Replace the coefficients of the CI states by the coefficients of the - eigenstates of the CI matrix - -`pt2_epstein_nesbet_sc2_projected `_ - compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution - .br - for the various N_st states, - .br - but with the correction in the denominator - .br - comming from the interaction of that determinant with all the others determinants - .br - that can be repeated by repeating all the double excitations - .br - : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) - .br - that could be repeated to this determinant. - .br - In addition, for the perturbative energetic contribution you have the standard second order - .br - e_2_pert = ^2/(Delta_E) - .br - and also the purely projected contribution - .br - H_pert_diag = c_pert - -`repeat_all_e_corr `_ - Undocumented - Needed Modules diff --git a/src/Selectors_full/README.rst b/src/Selectors_full/README.rst index 31bba35f..50010c0a 100644 --- a/src/Selectors_full/README.rst +++ b/src/Selectors_full/README.rst @@ -2,3 +2,96 @@ Selectors_full Module ===================== +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`coef_hf_selector `_ + energy of correlation per determinant respect to the Hartree Fock determinant + .br + for the all the double excitations in the selectors determinants + .br + E_corr_per_selectors(i) = * c(D_i)/c(HF) if |D_i> is a double excitation + .br + E_corr_per_selectors(i) = -1000.d0 if it is not a double excitation + .br + coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants + +`double_index_selectors `_ + degree of excitation respect to Hartree Fock for the wave function + .br + for the all the selectors determinants + .br + double_index_selectors = list of the index of the double excitations + .br + n_double_selectors = number of double excitations in the selectors determinants + +`e_corr_per_selectors `_ + energy of correlation per determinant respect to the Hartree Fock determinant + .br + for the all the double excitations in the selectors determinants + .br + E_corr_per_selectors(i) = * c(D_i)/c(HF) if |D_i> is a double excitation + .br + E_corr_per_selectors(i) = -1000.d0 if it is not a double excitation + .br + coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants + +`exc_degree_per_selectors `_ + degree of excitation respect to Hartree Fock for the wave function + .br + for the all the selectors determinants + .br + double_index_selectors = list of the index of the double excitations + .br + n_double_selectors = number of double excitations in the selectors determinants + +`n_double_selectors `_ + degree of excitation respect to Hartree Fock for the wave function + .br + for the all the selectors determinants + .br + double_index_selectors = list of the index of the double excitations + .br + n_double_selectors = number of double excitations in the selectors determinants + +`n_det_selectors `_ + For Single reference wave functions, the number of selectors is 1 : the + Hartree-Fock determinant + +`psi_selectors `_ + Determinants on which we apply for perturbation. + +`psi_selectors_coef `_ + Determinants on which we apply for perturbation. + +`psi_selectors_size `_ + Undocumented + +`threshold_selectors `_ + Percentage of the norm of the state-averaged wave function to + consider for the selectors + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Utils `_ + diff --git a/src/SingleRefMethod/README.rst b/src/SingleRefMethod/README.rst index c7100a95..837f5635 100644 --- a/src/SingleRefMethod/README.rst +++ b/src/SingleRefMethod/README.rst @@ -19,6 +19,9 @@ Documentation For Single reference wave functions, the generator is the Hartree-Fock determinant +`select_max `_ + Memo to skip useless selectors + Needed Modules diff --git a/src/SingleRefMethod/generators.irp.f b/src/SingleRefMethod/generators.irp.f index b12a681e..7cb1bc13 100644 --- a/src/SingleRefMethod/generators.irp.f +++ b/src/SingleRefMethod/generators.irp.f @@ -28,6 +28,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] call get_excitation_degree(HF_bitmask,psi_det(1,1,j),degree,N_int) if (degree == 0) then k = j + exit endif end do diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index b8a02470..cd71f1b7 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -48,7 +48,6 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) stop endif - !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(S_half,U,D,Vt,n,C,m) & !$OMP PRIVATE(i,j,k) @@ -178,7 +177,7 @@ subroutine apply_rotation(A,LDA,R,LDR,B,LDB,m,n) call dgemm('N','N',m,n,n,1.d0,A,LDA,R,LDR,0.d0,B,LDB) end -subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) +subroutine lapack_diagd(eigvalues,eigvectors,H,nmax,n) implicit none BEGIN_DOC ! Diagonalize matrix H @@ -244,3 +243,122 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) enddo deallocate(A,eigenvalues) end + +subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) + implicit none + BEGIN_DOC + ! Diagonalize matrix H + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax + double precision, intent(out) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) + double precision, intent(in) :: H(nmax,n) + double precision,allocatable :: eigenvalues(:) + double precision,allocatable :: work(:) + double precision,allocatable :: A(:,:) + integer :: lwork, info, i,j,l,k, liwork + + allocate(A(nmax,n),eigenvalues(n)) +! print*,'Diagonalization by jacobi' +! print*,'n = ',n + + A=H + lwork = 2*n*n + 6*n+ 1 + allocate (work(lwork)) + + lwork = -1 + call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + info ) + if (info < 0) then + print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int( work( 1 ) ) + deallocate (work) + + allocate (work(lwork)) + call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + info ) + deallocate(work) + + if (info < 0) then + print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'DSYEV Failed' + stop 1 + end if + + eigvectors = 0.d0 + eigvalues = 0.d0 + do j = 1, n + eigvalues(j) = eigenvalues(j) + do i = 1, n + eigvectors(i,j) = A(i,j) + enddo + enddo + deallocate(A,eigenvalues) +end + +subroutine lapack_partial_diag(eigvalues,eigvectors,H,nmax,n,n_st) + implicit none + BEGIN_DOC + ! Diagonalize matrix H + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax,n_st + double precision, intent(out) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) + double precision, intent(in) :: H(nmax,n) + double precision,allocatable :: work(:) + integer ,allocatable :: iwork(:), isuppz(:) + double precision,allocatable :: A(:,:) + integer :: lwork, info, i,j,l,k,m, liwork + + allocate(A(nmax,n)) + + A=H + lwork = 2*n*n + 6*n+ 1 + liwork = 5*n + 3 + allocate (work(lwork),iwork(liwork),isuppz(2*N_st)) + + lwork = -1 + liwork = -1 + call DSYEVR( 'V', 'I', 'U', n, A, nmax, 0.d0, 0.d0, 1, n_st, 1.d-10, m, eigvalues, eigvectors, nmax, isuppz, work, lwork, & + iwork, liwork, info ) + if (info < 0) then + print *, irp_here, ': DSYEVR: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int( work( 1 ) ) + liwork = iwork(1) + deallocate (work,iwork) + + allocate (work(lwork),iwork(liwork)) + call DSYEVR( 'V', 'I', 'U', n, A, nmax, 0.d0, 0.d0, 1, n_st, 1.d-10, m, eigvalues, eigvectors, nmax, isuppz, work, lwork, & + iwork, liwork, info ) + deallocate(work,iwork) + + if (info < 0) then + print *, irp_here, ': DSYEVR: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'DSYEVR Failed' + stop 1 + end if + + deallocate(A) +end diff --git a/src/Utils/README.rst b/src/Utils/README.rst index 3fd5bd8d..0ffe9840 100644 --- a/src/Utils/README.rst +++ b/src/Utils/README.rst @@ -92,7 +92,7 @@ Documentation into fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2) -`hermite `_ +`hermite `_ Hermite polynomial `multiply_poly `_ @@ -108,10 +108,10 @@ Documentation \int_0^1 dx \exp(-p x^2) x^n .br -`rint1 `_ +`rint1 `_ Standard version of rint -`rint_large_n `_ +`rint_large_n `_ Version of rint for large values of n `rint_sum `_