Merge branch 'dev' into csf

AUTHORS

@@ -0,0 +1,20 @@
+# If you contributed to this software, please make a pull request to add your
+# name to this list (alphabetical order of the last name)
+
+- Thomas Applencourt
+- Anouar Benali
+- Michel Caffarel
+- Grégoire David
+- Anthony Ferté
+- Yann Garniron
+- Kevin Gasperich
+- Vijay Gopal Chilkuri
+- Emmanuel Giner
+- Pierre-François Loos
+- Jean-Paul Malrieu
+- Julien Paquier
+- Barthélémy Pradines
+- Peter Reinhardt
+- Anthony Scemama
+- Julien Toulouse
+- Mickaël Véril

RELEASE_NOTES.org

@@ -0,0 +1,81 @@
+#+TITLE: Quantum Package Release notes
+
+* Version 2.2
+
+** New features
+
+** Changes
+
+  - Python3 replaces Python2
+  - Travis CI uses 3 jobs
+  - Moved Travis scripts into ~travis~ directory
+  - IRPF90 and EZFIO are now git submodules
+  - Now basis sets should be downloaded from basis-set-exchange website
+  - Added ~bse~ in the installable tools
+  - Documentation in ~src/README.rst~
+  - Added two-body reduced density matrix
+  - Added basis set correction
+  - Added CAS-based on-top density functional
+  - Improve PT2 computation for excited-states: Mostly 2x2
+    diagonalization, and some (n+1)x(n+1) diagonalizations
+  - Error bars for stochastic variance and norm of the perturbed wave function
+  - Improve PT2-matching for excited-states
+  - Compute the overlap of PT2 excited states
+  - Renamed SOP into CFG
+  - Improved parallelism in PT2 by splitting tasks
+  - Use max in multi-state PT2 instead of sum for the selection weight
+  - Added seniority
+  - Added excitation_max
+  - More tasks for distribueted Davidson
+  - Random guess vectors in Davidson have zeros to preserve symmetry
+  - Disk-based Davidson when too much memory is required
+  - Fixed bug in DIIS
+  - Fixed bug in molden (Au -> Angs)
+    
+*** User interface
+
+    - Added ~qp_basis~ script to install a basis set from the ~bse~
+      command-line tool
+    - Introduced ~n_det_qp_edit~, ~psi_det_qp_edit~, and
+ ~psi_coef_qp_edit~ to accelerate the opening of qp_edit with
+      large wave functions
+    - Removed ~etc/ninja.rc~
+    - Added flag to specify if the AOs are normalized 
+    - Added flag to specify if the primitive Gaussians are normalized
+    - Added ~lin_dep_cutoff~, the cutoff for linear dependencies
+    - Davidson convergence threshold can be adapted from PT2
+    - In ~density_for_dft~, ~no_core_density~ is now a logical
+    - Default for ~weight_selection~ has changed from 2 to 1
+    - Nullify_small_elements in matrices to keep symmetry
+    - Default of density functional changed from LDA to PBE
+    - Added ~no_vvvv_integrals~ flag
+    - Added ~pt2_min_parallel_tasks~ to control parallelism in PT2
+    - Added ~print_energy~
+    - Added ~print_hamiltonian~
+    - Added input for two body RDM
+
+*** Code
+
+    - Many bug fixes
+    - Changed electron-nucleus from ~e_n~ to ~n_e~ in names of variables
+    - Changed ~occ_pattern~ to ~configuration~
+    - Replaced ~List.map~ by a tail-recursive version ~Qputils.list_map~
+    - Added possible imaginary part in OCaml MO coefficients
+    - Added ~qp_clean_source_files.sh~ to remove non-ascii characters
+    - Added flag ~is_periodic~ for periodic systems
+    - Possibilities to handle complex integrals and complex MOs
+    - Moved pseuodpotential integrals out of ~ao_one_e_integrals~
+    - Removed Schwarz test and added logical functions
+ ~ao_two_e_integral_zero~ and ~ao_one_e_integral_zero~
+    - Introduced type for ~pt2_data~
+    - Banned excitations are used with far apart localized MOs
+    - S_z2_Sz is now included in S2
+    - S^2 in single precision
+    - Added Shank function
+    - Added utilities for periodic calculations
+
+    ao_one_e_integral_zero
+    banned_excitations
+      
+
+

TODO.org

@@ -1,29 +0,0 @@
-* popcnt pour avoir les determinants par NSOMO
-* Trier par Nsomo
-* Tableau (i) -> indice du 1er CFG qui a i SOMO
-* Boucler sur toutes les CFG mono-excitees par rapport a toutes les autres
-
-* p,q,cfg -> 0|1|2|3 (Type d'excitation)
-
-
-222200000000
-000010100101
-p->q
-222200200000
-000000000101
-
-SOMO->SOMO
-
-DOMO -> SOMO
-
-do p in DO
-  do q in not(DO or SO)
-    p->q + 2 DOMO->VMO
-
-
-do p in DO
-  do q in not(DO or SO)
-    p->q + 2 DOMO->VMO
-DOMO -> VMO
-test si bit=1 dans DO
-test si bit=1 dans VO

ocaml/qptypes_generator.ml

@@ -220,16 +220,16 @@ end = struct
   type t =
   | EN
   | HF
-  | SOP
+  | CFG
   [@@deriving sexp]
 
   let to_string = function
   | EN -> \"EN\"
   | HF -> \"HF\"
-  | SOP -> \"SOP\"
+  | CFG -> \"CFG\"
   let of_string  s =
     match (String.lowercase_ascii s) with
-    | \"sop\" -> SOP
+    | \"cfg\" -> CFG
     | \"en\"  -> EN
     | \"hf\"  -> HF
     | _ -> raise (Invalid_argument (\"Wrong Perturbation type : \"^s))

src/ao_two_e_ints/map_integrals.irp.f

@@ -698,81 +698,3 @@ subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values)
 end
 
 
-!subroutine dump_ao_integrals(filename)
-!  use map_module
-!  implicit none
-!  BEGIN_DOC
-!  ! Save to disk the |AO| integrals
-!  END_DOC
-!  character*(*), intent(in)      :: filename
-!  integer(cache_key_kind), pointer :: key(:)
-!  real(integral_kind), pointer   :: val(:)
-!  integer*8                      :: i,j, n
-!  if (.not.mpi_master) then
-!    return
-!  endif
-!  call ezfio_set_work_empty(.False.)
-!  open(unit=66,file=filename,FORM='unformatted')
-!  write(66) integral_kind, key_kind
-!  write(66) ao_integrals_map%sorted, ao_integrals_map%map_size,    &
-!      ao_integrals_map%n_elements
-!  do i=0_8,ao_integrals_map%map_size
-!    write(66) ao_integrals_map%map(i)%sorted, ao_integrals_map%map(i)%map_size,&
-!        ao_integrals_map%map(i)%n_elements
-!  enddo
-!  do i=0_8,ao_integrals_map%map_size
-!    key => ao_integrals_map%map(i)%key
-!    val => ao_integrals_map%map(i)%value
-!    n = ao_integrals_map%map(i)%n_elements
-!    write(66) (key(j), j=1,n), (val(j), j=1,n)
-!  enddo
-!  close(66)
-!
-!end
-
-
-!integer function load_ao_integrals(filename)
-!  implicit none
-!  BEGIN_DOC
-!  ! Read from disk the |AO| integrals
-!  END_DOC
-!  character*(*), intent(in)      :: filename
-!  integer*8                      :: i
-!  integer(cache_key_kind), pointer :: key(:)
-!  real(integral_kind), pointer   :: val(:)
-!  integer                        :: iknd, kknd
-!  integer*8                      :: n, j
-!  load_ao_integrals = 1
-!  open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
-!  read(66,err=98,end=98) iknd, kknd
-!  if (iknd /= integral_kind) then
-!    print *,  'Wrong integrals kind in file :', iknd
-!    stop 1
-!  endif
-!  if (kknd /= key_kind) then
-!    print *,  'Wrong key kind in file :', kknd
-!    stop 1
-!  endif
-!  read(66,err=98,end=98) ao_integrals_map%sorted, ao_integrals_map%map_size,&
-!      ao_integrals_map%n_elements
-!  do i=0_8, ao_integrals_map%map_size
-!    read(66,err=99,end=99) ao_integrals_map%map(i)%sorted,          &
-!        ao_integrals_map%map(i)%map_size, ao_integrals_map%map(i)%n_elements
-!    call cache_map_reallocate(ao_integrals_map%map(i),ao_integrals_map%map(i)%map_size)
-!  enddo
-!  do i=0_8, ao_integrals_map%map_size
-!    key => ao_integrals_map%map(i)%key
-!    val => ao_integrals_map%map(i)%value
-!    n = ao_integrals_map%map(i)%n_elements
-!    read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
-!  enddo
-!  call map_sort(ao_integrals_map)
-!  load_ao_integrals = 0
-!  return
-!  99 continue
-!  call map_deinit(ao_integrals_map)
-!  98 continue
-!  stop 'Problem reading ao_integrals_map file in work/'
-!
-!end
-!

src/ao_two_e_ints/two_e_integrals.irp.f

@@ -1166,13 +1166,13 @@ subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value)
   integer,intent(out)             :: n_integrals
   integer(key_kind),intent(out)   :: buffer_i(ao_num*ao_num)
   real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num)
+  logical, external               :: ao_two_e_integral_zero
 
   integer                         :: i,k
   double precision                :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2
   double precision                :: integral, wall_0
   double precision                :: thr
   integer                         :: kk, m, j1, i1
-  logical, external :: ao_two_e_integral_zero
 
   thr = ao_integrals_threshold
 

src/cipsi/selection.irp.f

@@ -541,6 +541,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   double precision, external :: diag_H_mat_elem_fock
   double precision :: E_shift
   double precision :: s_weight(N_states,N_states)
+  PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs
   do jstate=1,N_states
     do istate=1,N_states
       s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate))
@@ -592,6 +593,12 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       if(bannedOrb(p2, s2)) cycle
       if(banned(p1,p2)) cycle
 
+      if(pseudo_sym)then
+        if(dabs(mat(1, p1, p2)).lt.thresh_sym)then
+          w = 0.d0
+        endif
+      endif
+
       val = maxval(abs(mat(1:N_states, p1, p2)))
       if( val == 0d0) cycle
       call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
@@ -629,22 +636,35 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       endif
 
 
-      if (excitation_max >= 0) then
       integer :: degree
-        call get_excitation_degree(ref_bitmask(1,1),det(1,1),degree,N_int)
+      logical :: do_cycle
-        if (degree > excitation_max) cycle
+      if (excitation_max >= 0) then
+        do_cycle = .True.
+        do k=1,N_dominant_dets_of_cfgs
+          call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int)
+          do_cycle = do_cycle .and. (degree > excitation_max)
+        enddo
+        if (do_cycle) cycle
       endif
 
 
       if (excitation_alpha_max >= 0) then
-        call get_excitation_degree_spin(ref_bitmask(1,1),det(1,1),degree,N_int)
+        do_cycle = .True.
-        if (degree > excitation_alpha_max) cycle
+        do k=1,N_dominant_dets_of_cfgs
+          call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int)
+          do_cycle = do_cycle .and. (degree > excitation_alpha_max) 
+        enddo
+        if (do_cycle) cycle
       endif
 
 
       if (excitation_beta_max >= 0) then
-        call get_excitation_degree_spin(ref_bitmask(1,2),det(1,2),degree,N_int)
+        do_cycle = .True.
-        if (degree > excitation_beta_max) cycle
+        do k=1,N_dominant_dets_of_cfgs
+          call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int)
+          do_cycle = do_cycle .and. (degree > excitation_beta_max) 
+        enddo
+        if (do_cycle) cycle
       endif
 
 
@@ -780,12 +800,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
       end do
 
 
-      if(pseudo_sym)then
-        if(dabs(mat(1, p1, p2)).lt.thresh_sym)then
-          w = 0.d0
-        endif
-      endif
-
 !      w = dble(n) * w
 
 

src/determinants/configurations.irp.f

@@ -337,8 +337,6 @@ END_PROVIDER
 
 
 BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ]
- implicit none
- BEGIN_DOC
  ! Returns the index of the configuration for each determinant
  END_DOC
  integer :: i,j,k,r,l
@@ -542,3 +540,47 @@ end
 
 
 
+BEGIN_PROVIDER [ integer, dominant_cfg, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Configuration of the determinants with the largest weight, for each state
+ END_DOC
+ integer :: k
+ do k=1,N_states
+   dominant_cfg(k) = det_to_configuration(dominant_det(k))
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, N_dominant_dets_of_cfgs ]
+ implicit none
+ BEGIN_DOC
+ ! Number of determinants in all dominant determinants
+ END_DOC
+ integer                        :: k, sze
+
+ N_dominant_dets_of_cfgs = 0
+ do k=1,N_states
+   call configuration_to_dets_size(  &
+          psi_configuration(1,1,dominant_cfg(k)), &
+          sze, elec_alpha_num, N_int)
+   N_dominant_dets_of_cfgs += sze
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_dets_of_cfgs) ]
+ implicit none
+ BEGIN_DOC
+ ! Configuration of the determinants with the largest weight, for each state
+ END_DOC
+ integer :: i,k,sze
+ i=1
+ do k=1,N_states
+   sze = N_dominant_dets_of_cfgs
+   call configuration_to_dets( &
+          psi_configuration(1,1,dominant_cfg(k)), &
+          dominant_dets_of_cfgs(1,1,i), &
+          sze,elec_alpha_num,N_int)
+   i += sze
+ enddo
+END_PROVIDER

src/determinants/determinants.irp.f

@@ -256,6 +256,26 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ]
   enddo
 END_PROVIDER
 
+BEGIN_PROVIDER [ integer, dominant_det, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Determinant with the largest weight, for each state
+ END_DOC
+ integer :: i, k
+ double precision :: wmax, c
+ do k=1,N_states
+   wmax = 0.d0
+   do i=1,N_det
+     c = psi_coef(i,k)*psi_coef(i,k)
+     if (c > wmax) then
+       dominant_det(k) = i
+       wmax = c
+     endif
+   enddo
+ enddo
+
+END_PROVIDER
+
 
 
 !==============================================================================!