Merge pull request #11 from scemama/master

merge with scemama

plugins/MRCC_CASSD/EZFIO.cfg

@@ -7,11 +7,11 @@ interface: ezfio
 type: Threshold
 doc: Threshold on the convergence of the MRCC energy
 interface: ezfio,provider,ocaml
-default: 1.e-7
+default: 1.e-5
 
 [n_it_mrcc_max]
 type: Strictly_positive_int
 doc: Maximum number of MRCC iterations
 interface: ezfio,provider,ocaml
-default: 20
+default: 10
 

plugins/MRCC_CASSD/mrcc_cassd.irp.f

@@ -1,13 +1,46 @@
 program mrcc
   implicit none
-  if (.not.read_wf) then
+  read_wf = .True.
-    print *,  'read_wf has to be true.'
+  SOFT_TOUCH read_wf
-    stop 1
-  endif
   call print_cas_coefs
-  call run_mrcc
+  call set_generators_bitmasks_as_holes_and_particles
+  call run
 end
 
+subroutine run
+  implicit none
+  
+  integer :: i
+
+  double precision :: E_new, E_old, delta_e
+  integer :: iteration
+  double precision :: E_past(4), lambda
+  E_new = 0.d0
+  delta_E = 1.d0
+  iteration = 0
+  lambda = 1.d0
+  do while (delta_E > thresh_mrcc)
+    iteration += 1
+    print *,  '===========================' 
+    print *,  'MRCC Iteration', iteration
+    print *,  '===========================' 
+    print *,  ''
+    E_old = sum(ci_energy_dressed)
+    call write_double(6,ci_energy_dressed(1),"MRCC energy")
+    call diagonalize_ci_dressed(lambda)
+    E_new = sum(ci_energy_dressed)
+    delta_E = dabs(E_new - E_old)
+    call save_wavefunction
+    call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
+    if (iteration > n_it_mrcc_max) then
+      exit
+    endif
+  enddo
+  call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
+
+end
+
+
 subroutine print_cas_coefs
   implicit none
 
@@ -18,7 +51,7 @@ subroutine print_cas_coefs
     print *,  psi_cas_coef(i,:)
     call debug_det(psi_cas(1,1,i),N_int)
   enddo
-
   call write_double(6,ci_energy(1),"Initial CI energy")
+
 end
 

plugins/MRCC_CASSD/mrcc_noiter.irp.f

@@ -0,0 +1,41 @@
+program mrcc_noiter
+  implicit none
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call print_cas_coefs
+  call set_generators_bitmasks_as_holes_and_particles
+  call run
+end
+
+subroutine run
+  implicit none
+  
+  double precision :: lambda
+  integer :: i,j
+  do j=1,N_states_diag
+    do i=1,N_det
+      psi_coef(i,j) = CI_eigenvectors_dressed(i,j) 
+    enddo
+  enddo
+  SOFT_TOUCH psi_coef ci_energy_dressed
+  call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
+  call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
+  call save_wavefunction
+
+end
+
+
+subroutine print_cas_coefs
+  implicit none
+
+  integer :: i,j
+  print *,  'CAS'
+  print *,  '==='
+  do i=1,N_det_cas
+    print *,  psi_cas_coef(i,:)
+    call debug_det(psi_cas(1,1,i),N_int)
+  enddo
+  call write_double(6,ci_energy(1),"Initial CI energy")
+
+end
+

plugins/MRCC_Utils/mrcc_general.irp.f

@@ -11,12 +11,13 @@ subroutine mrcc_iterations
 
   double precision :: E_new, E_old, delta_e
   integer :: iteration,i_oscillations
-  double precision :: E_past(4)
+  double precision :: E_past(4), lambda
   E_new = 0.d0
   delta_E = 1.d0
   iteration = 0
   j = 1
   i_oscillations = 0
+  lambda = 1.d0
   do while (delta_E > 1.d-7)
     iteration += 1
     print *,  '===========================' 
@@ -25,10 +26,15 @@ subroutine mrcc_iterations
     print *,  ''
     E_old = sum(ci_energy_dressed)
     call write_double(6,ci_energy_dressed(1),"MRCC energy")
-    call diagonalize_ci_dressed
+    call diagonalize_ci_dressed(lambda)
     E_new = sum(ci_energy_dressed)
     delta_E = dabs(E_new - E_old)
-
+!    if (E_new > E_old) then
+!       lambda = lambda * 0.7d0
+!    else
+!       lambda = min(1.d0, lambda * 1.1d0)
+!    endif
+!    print *,  'energy lambda ', lambda
     E_past(j) = E_new
     j +=1
     call save_wavefunction

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -1,89 +1,45 @@
-! BEGIN_PROVIDER [ double precision, lambda_mrcc, (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),ihpsi_current(N_states)
-!  integer :: i_pert_count
-!  
-!  i_pert_count = 0
-!  lambda_mrcc = 0.d0
-!  
-!   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_current)
-!     do k=1,N_states
-!       if (ihpsi_current(k) == 0.d0) then
-!         ihpsi_current(k) = 1.d-32
-!       endif
-!       if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-6) then
-!         i_pert_count +=1
-!       else
-!         lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
-!       endif
-!     enddo
-!  enddo
-! 
-!  print*,'N_det_non_ref = ',N_det_non_ref
-!  print*,'Number of ignored determinants = ',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, (N_states,psi_det_size) ]
- implicit none
+  implicit none
- BEGIN_DOC
+  BEGIN_DOC
- ! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
+  ! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
- END_DOC
+  END_DOC
- integer :: i,k
+  integer :: i,k
   double precision               :: ihpsi_current(N_states)
   integer :: i_pert_count
   double precision :: hii, lambda_pert
-  
+
   i_pert_count = 0
   lambda_mrcc = 0.d0
-  
+
-   do i=1,N_det_non_ref
+    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, &
+      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)
+          size(psi_ref_coef,1), N_states,ihpsi_current)
-     call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
+      call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
-     do k=1,N_states
+      do k=1,N_states
-       if (ihpsi_current(k) == 0.d0) then
+        if (ihpsi_current(k) == 0.d0) then
-         ihpsi_current(k) = 1.d-32
+          ihpsi_current(k) = 1.d-32
-       endif
-       lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) 
-      if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-6 ) then
-           i_pert_count += 1
-              lambda_mrcc(k,i) = 0.d0
-!          lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
-!          if((ihpsi_current(k) * lambda_pert) < 0.5d0 * psi_non_ref_coef_restart(i,k) ) then
-!             lambda_mrcc(k,i) = 0.d0
-!          endif
-      endif
-        double precision, parameter :: x = 2.d0
-        if (lambda_mrcc(k,i) > x) then
-          lambda_mrcc(k,i) = x
-        else if (lambda_mrcc(k,i) < -x) then
-          lambda_mrcc(k,i) = -x
         endif
-     enddo
+        lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) )
-  enddo
+        lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
- 
+        if (lambda_pert / lambda_mrcc(k,i)  < 0.5d0 ) then
+          i_pert_count += 1
+          lambda_mrcc(k,i) = 0.d0
+        endif
+      enddo
+    enddo
+  
   print*,'N_det_non_ref = ',N_det_non_ref
-  print*,'Number of ignored determinants = ',i_pert_count
+  print*,'Number of ignored determinants = ',i_pert_count  
   print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
-  print*,'lambda min/max = ',maxval(dabs(lambda_mrcc)), minval(dabs(lambda_mrcc))
+  print*,'lambda max = ',maxval(dabs(lambda_mrcc))
 
 END_PROVIDER
 
 
-!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ]
+
-!implicit none
+
-!BEGIN_DOC
+
-!! Dressing matrix in SD basis
+
-!END_DOC
-!delta_ij_non_ref = 0.d0
-!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref)
-!END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
  implicit none
@@ -218,17 +174,19 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
 
 END_PROVIDER
 
-subroutine diagonalize_CI_dressed
+subroutine diagonalize_CI_dressed(lambda)
   implicit none
   BEGIN_DOC
 !  Replace the coefficients of the CI states by the coefficients of the 
 !  eigenstates of the CI matrix
   END_DOC
+  double precision, intent(in) :: lambda
   integer :: i,j
   do j=1,N_states_diag
     do i=1,N_det
-      psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
+      psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j) 
     enddo
+    call normalize(psi_coef(1,j), N_det)
   enddo
   SOFT_TOUCH psi_coef 
 

plugins/Psiref_CAS/psi_ref.irp.f

@@ -23,6 +23,7 @@ use bitmasks
       psi_ref_coef(i,k) = psi_cas_coef(i,k)
     enddo
   enddo
+  call normalize(psi_ref_coef, N_det_ref)
 
 END_PROVIDER
 

tests/bats/qp.bats

@@ -53,10 +53,10 @@ function test_exe() {
 }
 
 function run_HF() {
-  thresh=1.e-8
+  thresh=1.e-7
   test_exe SCF || skip
   ezfio set_file $1
-  ezfio set hartree_fock thresh_scf 1.e-10
+  ezfio set hartree_fock thresh_scf 1.e-11
   qp_run SCF $1 
   energy="$(ezfio get hartree_fock energy)"
   eq $energy $2 $thresh
@@ -155,7 +155,7 @@ function run_all_1h_1p() {
   ezfio set determinants read_wf True
   qp_run mrcc_cassd $INPUT  
   energy="$(ezfio get mrcc_cassd energy)"
-  eq $energy -76.2289109271715     1.E-3
+  eq $energy -76.2286679037553  1.e-4
   
 }
 
@@ -166,7 +166,7 @@ function run_all_1h_1p() {
 }
 
 @test "SCF H2O VDZ pseudo" {
-  run_HF  h2o_pseudo.ezfio  -0.169483703904991E+02
+  run_HF  h2o_pseudo.ezfio  -16.9483708495521
 }
 
 @test "FCI H2O VDZ pseudo" {