MRCC Manu

ocaml/qp_create_ezfio_from_xyz.ml

@@ -14,7 +14,7 @@ let spec =
   +> flag "m" (optional_with_default 1 int)
      ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1."
   +> flag "p" no_arg
-     ~doc:"Using pseudopotentials"
+     ~doc:" Using pseudopotentials"
   +> anon ("xyz_file" %: string)
 ;;
 

plugins/MRCC_CASSD/.gitignore

@@ -14,7 +14,7 @@ Integrals_Bielec
 Integrals_Monoelec
 MOGuess
 MO_Basis
-MRCC_Utils_new
+MRCC_Utils
 Makefile
 Makefile.depend
 Nuclei

plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils_new
+Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils

plugins/MRCC_CASSD/README.rst

@@ -16,7 +16,7 @@ Needed Modules
 * `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
 * `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
 * `Psiref_CAS <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_CAS>`_
-* `MRCC_Utils_new <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils_new>`_
+* `MRCC_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils>`_
 
 Documentation
 =============

plugins/MRCC_Utils/README.rst

@@ -21,19 +21,19 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. by the `update_README.py` script.
 
-`ci_eigenvectors_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L76>`_
+`ci_eigenvectors_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L166>`_
   Eigenvectors/values of the CI matrix
 
 
-`ci_eigenvectors_s2_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L77>`_
+`ci_eigenvectors_s2_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L167>`_
   Eigenvectors/values of the CI matrix
 
 
-`ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L75>`_
+`ci_electronic_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L165>`_
   Eigenvectors/values of the CI matrix
 
 
-`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L142>`_
+`ci_energy_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L232>`_
   N_states lowest eigenvalues of the dressed CI matrix
 
 
@@ -77,15 +77,15 @@ Documentation
   Initial guess vectors are not necessarily orthonormal
 
 
-`delta_ii <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L39>`_
+`delta_ii <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L104>`_
   Dressing matrix in N_det basis
 
 
-`delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L38>`_
+`delta_ij <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L103>`_
   Dressing matrix in N_det basis
 
 
-`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L157>`_
+`diagonalize_ci_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L247>`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
@@ -111,7 +111,7 @@ Documentation
   Assume N_int is already provided.
 
 
-`h_matrix_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L50>`_
+`h_matrix_dressed <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L140>`_
   Dressed H with Delta_ij
 
 
@@ -123,11 +123,15 @@ Documentation
   H_jj : array of <j|H|j>
 
 
-`lambda_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L1>`_
+`lambda_mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L5>`_
   cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
 
 
-`lambda_pert <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L2>`_
+`lambda_mrcc_tmp <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L81>`_
+  Undocumented
+
+
+`lambda_pert <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L6>`_
   cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
 
 
@@ -143,6 +147,14 @@ Documentation
   Undocumented
 
 
+`oscillations <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L86>`_
+  Undocumented
+
+
+`pert_determinants <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_utils.irp.f#L1>`_
+  Undocumented
+
+
 `psi_ref_lock <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_dress.irp.f#L3>`_
   Locks on ref determinants to fill delta_ij
 
@@ -151,6 +163,6 @@ Documentation
   Undocumented
 
 
-`set_generators_bitmasks_as_holes_and_particles <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L39>`_
+`set_generators_bitmasks_as_holes_and_particles <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils/mrcc_general.irp.f#L69>`_
   Undocumented
 

plugins/MRCC_Utils/mrcc_general.irp.f

@@ -10,11 +10,14 @@ subroutine mrcc_iterations
   integer :: i,j
 
   double precision :: E_new, E_old, delta_e
-  integer :: iteration
+  integer :: iteration,i_oscillations
+  double precision :: E_past(4)
   E_new = 0.d0
   delta_E = 1.d0
   iteration = 0
-  do while (delta_E > 100.d0*davidson_threshold)
+  j = 1
+  i_oscillations = 0
+  do while (delta_E > 1.d-7)
     iteration += 1
     print *,  '===========================' 
     print *,  'MRCC Iteration', iteration
@@ -25,10 +28,37 @@ subroutine mrcc_iterations
     call diagonalize_ci_dressed
     E_new = sum(ci_energy_dressed)
     delta_E = dabs(E_new - E_old)
-    if (iteration > 20) then
-      exit
+
+    E_past(j) = E_new
+    j +=1
+    if(j>4)then
+     j=1
+    endif
+    if(iteration > 4) then 
+     if(delta_E > 1.d-10)then
+      if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then 
+       print*,'OSCILLATIONS !!!'
+       oscillations = .True.
+       i_oscillations +=1
+       lambda_mrcc_tmp = lambda_mrcc
+      endif
+     endif
     endif
     call save_wavefunction
+!   if (i_oscillations > 5) then
+!    exit
+!   endif
+    if (iteration > 200) then
+      exit
+    endif
+    print*,'------------'
+    print*,'VECTOR'
+    do i = 1, N_det_ref
+     print*,''
+     print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1)
+     print*,'delta_ii(i,1)     = ',delta_ii(i,1)
+    enddo
+    print*,'------------'
   enddo
   call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
   call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -1,26 +1,91 @@
+ BEGIN_PROVIDER [integer, pert_determinants, (N_states, psi_det_size) ]
+ END_PROVIDER 
+
+
  BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
 &BEGIN_PROVIDER [ double precision, lambda_pert, (N_states,psi_det_size) ] 
  implicit none
  BEGIN_DOC
  ! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
  END_DOC
- integer :: i,k
- double precision               :: ihpsi(N_states), hii
+ integer :: i,k,j
+ double precision               :: ihpsi(N_states), hii,delta_e_eff,ihpsi_current(N_states),hij
+ integer :: i_ok,i_pert,i_pert_count
+ i_ok = 0
+
+ double precision :: phase_restart(N_states),tmp
+ do k = 1, N_states
+  phase_restart(k) = dsign(1.d0,psi_ref_coef_restart(1,k)/psi_ref_coef(1,k))
+ enddo
+ i_pert_count = 0
  
  do i=1,N_det_non_ref
-   call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
-       size(psi_ref_coef,1), n_states, ihpsi)
-   call i_h_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
+   call i_h_psi(psi_non_ref(1,1,i), psi_ref_restart, psi_ref_coef_restart, N_int, N_det_ref,&
+       size(psi_ref_coef_restart,1), n_states, ihpsi)
+   call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
    do k=1,N_states
      lambda_pert(k,i) = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
-     if (dabs(ihpsi(k)).le.1.d-3) then
+     call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
+     tmp = psi_non_ref_coef(i,k)/ihpsi_current(k)
+     i_pert = 1
+     if((ihpsi(k) * lambda_pert(k,i))/psi_non_ref_coef_restart(i,k) .ge. 0.5d0 & 
+        .and. (ihpsi(k) * lambda_pert(k,i))/psi_non_ref_coef_restart(i,k) > 0.d0 )then  ! test on the first order coefficient
+      i_pert = 0
+     endif
+     do j = 1, N_det_ref
+      call i_H_j(psi_non_ref(1,1,i),psi_ref(1,1,j),N_int,hij)
+      if(dabs(hij * tmp).ge.0.5d0)then
+       i_pert_count +=1
+       i_pert = 1
+       exit
+      endif
+     enddo
+     if( i_pert == 1)then
+      pert_determinants(k,i) = i_pert
+     endif
+     if(pert_determinants(k,i) == 1)then
+       i_ok +=1
        lambda_mrcc(k,i) = lambda_pert(k,i)
      else
-       lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi(k)
+       lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
      endif
    enddo
  enddo
+!if(oscillations)then
+! print*,'AVERAGING the lambda_mrcc with those of the previous iterations'
+! do i = 1, N_det_non_ref
+!  do k = 1, N_states
 
+!   double precision :: tmp
+!   tmp = lambda_mrcc(k,i)
+!   lambda_mrcc(k,i) += lambda_mrcc_tmp(k,i)
+!   lambda_mrcc(k,i) = lambda_mrcc(k,i) * 0.5d0
+!   if(dabs(tmp - lambda_mrcc(k,i)).ge.1.d-9)then
+!   print*,''
+!   print*,'i = ',i
+!   print*,'psi_non_ref_coef(i,k) = ',psi_non_ref_coef(i,k)
+!   print*,'lambda_mrcc(k,i)     = ',lambda_mrcc(k,i)
+!   print*,'                 tmp = ',tmp
+!   endif
+!  enddo
+! enddo
+!endif
+ print*,'N_det_non_ref = ',N_det_non_ref
+ print*,'Number of Perturbatively treated determinants = ',i_ok
+ print*,'i_pert_count = ',i_pert_count
+ print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, lambda_mrcc_tmp, (N_states,psi_det_size) ]
+ implicit none
+ lambda_mrcc_tmp = 0.d0
+END_PROVIDER 
+
+BEGIN_PROVIDER [ logical, oscillations ]
+ implicit none
+ oscillations = .False.
 END_PROVIDER 
 
 
@@ -45,6 +110,31 @@ END_PROVIDER
  delta_ij = 0.d0
  delta_ii = 0.d0
  call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref)
+ double precision :: max_delta
+ double precision :: accu
+ integer :: imax,jmax
+ max_delta = 0.d0
+ accu = 0.d0
+ do i = 1, N_det_ref
+  do j = 1, N_det_non_ref
+   accu += psi_non_ref_coef(j,1) * psi_ref_coef(i,1) * delta_ij(i,j,1)
+   if(dabs(delta_ij(i,j,1)).gt.max_delta)then
+    max_delta = dabs(delta_ij(i,j,1))
+    imax = i
+    jmax = j
+   endif
+  enddo
+ enddo
+ print*,''
+ print*,''
+ print*,'<psi| Delta H |psi> = ',accu
+ print*,'MAX VAL OF DRESING = ',delta_ij(imax,jmax,1)
+ print*,'imax,jmax = ',imax,jmax
+ print*,'psi_ref_coef(imax,1)     = ',psi_ref_coef(imax,1)
+ print*,'psi_non_ref_coef(jmax,1) = ',psi_non_ref_coef(jmax,1)
+ do i = 1, N_det_ref
+  print*,'delta_ii(i,1)     = ',delta_ii(i,1)
+ enddo
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]

plugins/Psiref_CAS/README.rst

@@ -31,6 +31,14 @@ Documentation
   CAS wave function, defined from the application of the CAS bitmask on the
   determinants. idx_cas gives the indice of the CAS determinant in psi_det.
 
+
+`psi_ref_coef_restart <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_CAS/psi_ref.irp.f#L30>`_
+  Projection of the CAS wave function on the restart wave function.
+
+
+`psi_ref_restart <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_CAS/psi_ref.irp.f#L29>`_
+  Projection of the CAS wave function on the restart wave function.
+
 Needed Modules
 ==============
 

plugins/Psiref_CAS/psi_ref.irp.f

@@ -26,3 +26,29 @@ use bitmasks
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [ integer(bit_kind), psi_ref_restart, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_ref_coef_restart,  (psi_det_size,n_states) ]
+  implicit none
+  BEGIN_DOC
+  ! Projection of the CAS wave function on the restart wave function. 
+  END_DOC
+  integer :: i,j,k
+  integer, save                  :: ifirst
+
+  if(ifirst == 0)then
+   ifirst = 1
+   do i=1,N_det_ref
+     do k=1,N_int
+       psi_ref_restart(k,1,i) = psi_cas(k,1,i)
+       psi_ref_restart(k,2,i) = psi_cas(k,2,i)
+     enddo
+   enddo
+   do k=1,N_states
+     do i=1,N_det_ref
+       psi_ref_coef_restart(i,k) = psi_cas_coef(i,k)
+     enddo
+   enddo
+  endif
+
+END_PROVIDER
+

plugins/Psiref_Utils/README.rst

@@ -13,11 +13,11 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. by the `update_README.py` script.
 
-`get_index_in_psi_ref_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L140>`_
+`get_index_in_psi_ref_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L182>`_
   Returns the index of the determinant in the ``psi_ref_sorted_bit`` array
 
 
-`h_matrix_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L74>`_
+`h_matrix_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L116>`_
   Undocumented
 
 
@@ -42,7 +42,7 @@ Documentation
   idx_non_ref_rev gives the reverse.
 
 
-`is_in_psi_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L126>`_
+`is_in_psi_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L168>`_
   True if the determinant ``det`` is in the wave function
 
 
@@ -60,7 +60,7 @@ Documentation
   been done going from psi_ref to psi_non_ref
 
 
-`psi_coef_ref_diagonalized <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L86>`_
+`psi_coef_ref_diagonalized <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L128>`_
   Undocumented
 
 
@@ -78,12 +78,26 @@ Documentation
   idx_non_ref_rev gives the reverse.
 
 
-`psi_non_ref_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L62>`_
+`psi_non_ref_coef_restart <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L62>`_
+  Set of determinants which are not part of the reference, defined from the application
+  of the reference bitmask on the determinants.
+  idx_non_ref gives the indice of the determinant in psi_det.
+  But this is with respect to the restart wave function.
+
+
+`psi_non_ref_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L104>`_
   Reference determinants sorted to accelerate the search of a random determinant in the wave
   function.
 
 
-`psi_non_ref_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L61>`_
+`psi_non_ref_restart <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L61>`_
+  Set of determinants which are not part of the reference, defined from the application
+  of the reference bitmask on the determinants.
+  idx_non_ref gives the indice of the determinant in psi_det.
+  But this is with respect to the restart wave function.
+
+
+`psi_non_ref_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L103>`_
   Reference determinants sorted to accelerate the search of a random determinant in the wave
   function.
 
@@ -93,11 +107,11 @@ Documentation
   function.
 
 
-`psi_ref_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L105>`_
+`psi_ref_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L147>`_
   Undocumented
 
 
-`psi_ref_energy_diagonalized <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L87>`_
+`psi_ref_energy_diagonalized <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_Utils/psi_ref_utils.irp.f#L129>`_
   Undocumented
 
 

plugins/Psiref_Utils/psi_ref_utils.irp.f

@@ -58,6 +58,48 @@ END_PROVIDER
  N_det_non_ref = i_non_ref
 END_PROVIDER
 
+ BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref_restart,  (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_non_ref_coef_restart, (psi_det_size,n_states) ]
+ implicit none
+ BEGIN_DOC
+  ! Set of determinants which are not part of the reference, defined from the application
+  ! of the reference bitmask on the determinants. 
+  ! idx_non_ref gives the indice of the determinant in psi_det.
+  ! But this is with respect to the restart wave function. 
+ END_DOC
+ integer                        :: i_non_ref,j,k
+ integer                        :: degree
+ logical                        :: in_ref
+ integer, save                  :: ifirst = 0
+ if(ifirst==0)then
+  ifirst = 1
+  i_non_ref =0
+  do k=1,N_det
+    in_ref = .False.
+    do j=1,N_det_ref
+      call get_excitation_degree(psi_ref(1,1,j), psi_det(1,1,k), degree, N_int)
+      if (degree == 0) then
+        in_ref = .True.
+        exit
+      endif
+    enddo
+    if (.not.in_ref) then
+      double precision :: hij
+      i_non_ref += 1
+      do j=1,N_int
+        psi_non_ref_restart(j,1,i_non_ref) = psi_det(j,1,k)
+        psi_non_ref_restart(j,2,i_non_ref) = psi_det(j,2,k)
+      enddo
+      do j=1,N_states
+        psi_non_ref_coef_restart(i_non_ref,j) = psi_coef(k,j)
+      enddo
+    endif
+  enddo
+ endif
+END_PROVIDER
+
+
+
  BEGIN_PROVIDER [ integer(bit_kind), psi_non_ref_sorted_bit, (N_int,2,psi_det_size) ]
 &BEGIN_PROVIDER [ double precision, psi_non_ref_coef_sorted_bit, (psi_det_size,N_states) ]
  implicit none