Merge branch 'develop'

ocaml/qp_create_guess.ml

@@ -88,8 +88,9 @@ let run ~multiplicity ezfio_file =
       ~alpha:(Elec_alpha_number.of_int alpha_new)
       ~beta:(Elec_beta_number.of_int beta_new) pair )
   in
+
   let c = 
-    Array.create ~len:(List.length determinants) (Det_coef.of_float 1.)
+    Array.init (List.length determinants) (fun _ -> Det_coef.of_float ((Random.float 2.)-.1.))
   in
  
   determinants 

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -751,6 +751,10 @@ END_PROVIDER
     end do
     deallocate(lref)
 
+    do i=1,N_det_non_ref
+      rho_mrcc(i,s) = rho_mrcc_init(i) 
+    enddo
+
     x_new = x
     
     double precision               :: factor, resold
@@ -758,14 +762,8 @@ END_PROVIDER
     resold = huge(1.d0)
 
     do k=0,10*hh_nex
-      !$OMP PARALLEL default(shared) private(cx, i, a_col, a_coll)
+      res = 0.d0
-      
+      !$OMP PARALLEL default(shared) private(cx, i, a_col, a_coll) reduction(+:res)
-      !$OMP DO
-      do i=1,N_det_non_ref
-        rho_mrcc(i,s) = rho_mrcc_init(i) 
-      enddo
-      !$OMP END DO
-      
       !$OMP DO
       do a_coll = 1, n_exc_active
         a_col = active_pp_idx(a_coll)
@@ -774,23 +772,12 @@ END_PROVIDER
           cx = cx + x(mrcc_AtA_ind(i)) * mrcc_AtA_val(s,i)
         end do
         x_new(a_col) = AtB(a_col) + cx * factor
-      end do
-      !$OMP END DO
-
-      !$OMP END PARALLEL
-      
-      
-      res = 0.d0
-      do a_coll=1,n_exc_active
-        a_col = active_pp_idx(a_coll)
-        do j=1,N_det_non_ref
-          i = active_excitation_to_determinants_idx(j,a_coll)
-          if (i==0) exit
-          rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * X_new(a_col)
-        enddo
         res = res + (X_new(a_col) - X(a_col))*(X_new(a_col) - X(a_col))
         X(a_col) = X_new(a_col)
       end do
+      !$OMP END DO
+      !$OMP END PARALLEL
+      
       if (res > resold) then
         factor = factor * 0.5d0
       endif
@@ -802,6 +789,42 @@ END_PROVIDER
       
       if(res < 1d-10) exit
     end do
+    dIj_unique(1:size(X), s) = X(1:size(X))
+
+!    double precision, external :: ddot
+!    if (ddot (size(X), dIj_unique, 1, X, 1) < 0.d0) then
+!      dIj_unique(1:size(X),s) = -X(1:size(X))
+!    endif
+    
+  enddo
+
+  ! Adjust phase of dIj_unique
+
+!  double precision :: snorm
+!  X = 0.d0
+!  snorm = 0.d0
+!  do s=1,N_states
+!    norm = 0.d0
+!    do i=1,N_det_non_ref
+!      norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s)
+!    enddo
+!    norm = dsqrt(norm)
+!    X(1:size(X)) = X(1:size(X)) + dIj_unique(1:size(X),s) * norm
+!    snorm += norm
+!  enddo
+!  X = X/snorm
+
+  do s=1,N_states
+
+    do a_coll=1,n_exc_active
+      a_col = active_pp_idx(a_coll)
+      do j=1,N_det_non_ref
+        i = active_excitation_to_determinants_idx(j,a_coll)
+        if (i==0) exit
+        rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * dIj_unique(a_col,s)
+!        rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * X(a_col)
+      enddo
+    end do
 
     norm = 0.d0
     do i=1,N_det_non_ref
@@ -814,122 +837,11 @@ END_PROVIDER
     enddo
     ! Norm now contains the norm of Psi + A.X
     
-    print *, k, "res : ", res, "norm : ", sqrt(norm)
+    print *, "norm : ", sqrt(norm)
+   enddo
      
-!---------------
-! double precision               :: e_0, overlap
-! double precision, allocatable  :: u_0(:)
-! integer(bit_kind), allocatable :: keys_tmp(:,:,:)
-! allocate (u_0(N_det), keys_tmp(N_int,2,N_det) )
-! k=0
-! overlap = 0.d0
-! do i=1,N_det_ref
-!   k = k+1
-!   u_0(k) = psi_ref_coef(i,1)
-!   keys_tmp(:,:,k) = psi_ref(:,:,i)
-!   overlap += u_0(k)*psi_ref_coef(i,1)
-! enddo
-! norm = 0.d0
-! do i=1,N_det_non_ref
-!   k = k+1
-!   u_0(k) = psi_non_ref_coef(i,1)
-!   keys_tmp(:,:,k) = psi_non_ref(:,:,i)
-!   overlap += u_0(k)*psi_non_ref_coef(i,1)
-! enddo
-! 
-! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
-! print *,  'Energy of |Psi_CASSD> : ', e_0 + nuclear_repulsion, overlap
-!
-! k=0
-! overlap = 0.d0
-! do i=1,N_det_ref
-!   k = k+1
-!   u_0(k) = psi_ref_coef(i,1)
-!   keys_tmp(:,:,k) = psi_ref(:,:,i)
-!   overlap += u_0(k)*psi_ref_coef(i,1)
-! enddo
-! norm = 0.d0
-! do i=1,N_det_non_ref
-!   k = k+1
-!   ! f is such that f.\tilde{c_i} = c_i
-!   f = psi_non_ref_coef(i,1) / rho_mrcc(i,1)
-!   
-!   ! Avoid numerical instabilities
-!   f = min(f,2.d0)
-!   f = max(f,-2.d0)
-!
-!   f = 1.d0
-!
-!   u_0(k) = rho_mrcc(i,1)*f
-!   keys_tmp(:,:,k) = psi_non_ref(:,:,i)
-!   norm += u_0(k)**2
-!   overlap += u_0(k)*psi_non_ref_coef(i,1)
-! enddo
-! 
-! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
-! print *,  'Energy of |(1+T)Psi_0> : ', e_0 + nuclear_repulsion, overlap
-!
-! f = 1.d0/norm
-! norm = 1.d0
-! do i=1,N_det_ref
-!  norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s)
-! enddo
-! f = dsqrt(f*norm)
-! overlap = norm
-! do i=1,N_det_non_ref
-!   u_0(k) = rho_mrcc(i,1)*f
-!   overlap += u_0(k)*psi_non_ref_coef(i,1)
-! enddo
-!
-! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
-! print *,  'Energy of |(1+T)Psi_0> (normalized) : ', e_0 + nuclear_repulsion,  overlap
-!
-! k=0
-! overlap = 0.d0
-! do i=1,N_det_ref
-!   k = k+1
-!   u_0(k) = psi_ref_coef(i,1)
-!   keys_tmp(:,:,k) = psi_ref(:,:,i)
-!   overlap += u_0(k)*psi_ref_coef(i,1)
-! enddo
-! norm = 0.d0
-! do i=1,N_det_non_ref
-!   k = k+1
-!   ! f is such that f.\tilde{c_i} = c_i
-!   f = psi_non_ref_coef(i,1) / rho_mrcc(i,1)
-!   
-!   ! Avoid numerical instabilities
-!   f = min(f,2.d0)
-!   f = max(f,-2.d0)
-!
-!   u_0(k) = rho_mrcc(i,1)*f
-!   keys_tmp(:,:,k) = psi_non_ref(:,:,i)
-!   norm += u_0(k)**2
-!   overlap += u_0(k)*psi_non_ref_coef(i,1)
-! enddo
-! 
-! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
-! print *,  'Energy of |(1+T)Psi_0> (mu_i): ', e_0 + nuclear_repulsion, overlap
-!
-! f = 1.d0/norm
-! norm = 1.d0
-! do i=1,N_det_ref
-!  norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s)
-! enddo
-! overlap = norm
-! f = dsqrt(f*norm)
-! do i=1,N_det_non_ref
-!   u_0(k) = rho_mrcc(i,1)*f
-!   overlap += u_0(k)*psi_non_ref_coef(i,1)
-! enddo
-!
-! call u_0_H_u_0(e_0,u_0,N_det,keys_tmp,N_int,1,N_det)
-! print *,  'Energy of |(1+T)Psi_0> (normalized mu_i) : ', e_0 + nuclear_repulsion, overlap
-!
-! deallocate(u_0, keys_tmp)
-!
-!---------------
         
+   do s=1,N_states
      norm = 0.d0
      double precision               :: f
      do i=1,N_det_non_ref
@@ -937,12 +849,16 @@ END_PROVIDER
          rho_mrcc(i,s) = 1.d-32
        endif
 
+       if (lambda_type == 2) then
+         f = 1.d0
+       else
         ! f is such that f.\tilde{c_i} = c_i
         f = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
 
         ! Avoid numerical instabilities
         f = min(f,2.d0)
         f = max(f,-2.d0)
+      endif
 
        norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s)
        rho_mrcc(i,s) = f
@@ -977,7 +893,6 @@ END_PROVIDER
      ! rho_mrcc now contains the product of the scaling factors and the
      ! normalization constant
     
-    dIj_unique(1:size(X), s) = X(1:size(X))
   end do
 
 END_PROVIDER
@@ -1018,6 +933,7 @@ double precision function get_dij_index(II, i, s, Nint)
   else if(lambda_type == 2) then
     call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int)
     get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase
+    get_dij_index = get_dij_index * rho_mrcc(i,s) 
   end if
 end function