Removed eI

src/cipsi/pt2_stoch_routines.irp.f

@@ -242,8 +242,8 @@ subroutine ZMQ_pt2(E, pt2_data, relative_error, N_in)
       mem_collector = 8.d0 *                  & ! bytes
             ( 1.d0*pt2_n_tasks_max            & ! task_id, index
             + 0.635d0*N_det_generators        & ! f,d
-            + 3.d0*N_det_generators*N_states  & ! eI, vI, nI
             + pt2_n_tasks_max*pt2_type_size(N_states)/8 & ! pt2_data_task
+            + N_det_generators*pt2_type_size(N_states)/8  & ! pt2_data_I
             + 4.d0*(pt2_N_teeth+1)            & ! S, S2, T2, T3
             + 1.d0*(N_int*2.d0*N + N)         & ! selection buffer
             + 1.d0*(N_int*2.d0*N + N)         & ! sort selection buffer
@@ -360,9 +360,10 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
   integer, intent(in)            :: N_
 
   type(pt2_type), allocatable    :: pt2_data_task(:)
-  double precision, allocatable      :: eI(:,:), Se(:), Se2(:)
+  type(pt2_type), allocatable    :: pt2_data_I(:)
-  double precision, allocatable      :: vI(:,:), Sv(:), Sv2(:)
+  type(pt2_type), allocatable    :: pt2_data_S(:)
-  double precision, allocatable      :: nI(:,:), Sn(:), Sn2(:)
+  type(pt2_type), allocatable    :: pt2_data_S2(:)
+  type(pt2_type)                 :: pt2_data_teeth
   integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket
   integer(ZMQ_PTR)               :: zmq_to_qp_run_socket
   integer, external :: zmq_delete_tasks_async_send
@@ -370,6 +371,8 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
   integer, external :: zmq_abort
   integer, external :: pt2_find_sample_lr
 
+  PROVIDE pt2_stoch_istate
+
   integer :: more, n, i, p, c, t, n_tasks, U
   integer, allocatable :: task_id(:)
   integer, allocatable :: index(:)
@@ -400,12 +403,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
   allocate(task_id(pt2_n_tasks_max), index(pt2_n_tasks_max), f(N_det_generators))
   allocate(d(N_det_generators+1))
   allocate(pt2_data_task(pt2_n_tasks_max))
-  allocate(eI(N_states, N_det_generators))
+  allocate(pt2_data_I(N_det_generators))
-  allocate(vI(N_states, N_det_generators))
+  allocate(pt2_data_S(pt2_N_teeth+1))
-  allocate(nI(N_states, N_det_generators))
+  allocate(pt2_data_S2(pt2_N_teeth+1))
-  allocate(Se(pt2_N_teeth+1), Se2(pt2_N_teeth+1))
-  allocate(Sv(pt2_N_teeth+1), Sv2(pt2_N_teeth+1))
-  allocate(Sn(pt2_N_teeth+1), Sn2(pt2_N_teeth+1))
 
 
 
@@ -417,18 +417,19 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
   pt2_data % pt2_err(pt2_stoch_istate) = huge(1.)
   pt2_data % variance(pt2_stoch_istate) = huge(1.)
   pt2_data % norm2(pt2_stoch_istate) = 0.d0
-  Se(:) = 0d0
-  Sv(:) = 0d0
-  Sn(:) = 0d0
-  Se2(:) = 0d0
-  Sv2(:) = 0d0
-  Sn2(:) = 0d0
   n = 1
   t = 0
   U = 0
-  eI(:,:) = 0d0
+  do i=1,pt2_n_tasks_max
-  vI(:,:) = 0d0
+    call pt2_alloc(pt2_data_task(i),N_states)
-  nI(:,:) = 0d0
+  enddo
+  do i=1,pt2_N_teeth+1
+    call pt2_alloc(pt2_data_S(i),N_states)
+    call pt2_alloc(pt2_data_S2(i),N_states)
+  enddo
+  do i=1,N_det_generators
+    call pt2_alloc(pt2_data_I(i),N_states)
+  enddo
   f(:) = pt2_F(:)
   d(:) = .false.
   n_tasks = 0
@@ -456,9 +457,9 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
           v0 = 0.d0
           n0 = 0.d0
           do i=pt2_n_0(t),1,-1
-            E0 += eI(pt2_stoch_istate, i)
+            E0 += pt2_data_I(i) % pt2(pt2_stoch_istate)
-            v0 += vI(pt2_stoch_istate, i)
+            v0 += pt2_data_I(i) % variance(pt2_stoch_istate)
-            n0 += nI(pt2_stoch_istate, i)
+            n0 += pt2_data_I(i) % norm2(pt2_stoch_istate)
           end do
         else
           exit
@@ -468,26 +469,19 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
       ! Add Stochastic part
       c = pt2_R(n)
       if(c > 0) then
-!print *,  'c>0'
+
-        x  = 0d0
+        call pt2_alloc(pt2_data_teeth,N_states)
-        x2 = 0d0
-        x3 = 0d0
         do p=pt2_N_teeth, 1, -1
           v = pt2_u_0 + pt2_W_T * (pt2_u(c) + dble(p-1))
           i = pt2_find_sample_lr(v, pt2_cW,pt2_n_0(p),pt2_n_0(p+1))
-          x  += eI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i)
+          call pt2_add ( pt2_data_teeth,  pt2_W_T / pt2_w(i),  pt2_data_I(i) )
-          x2 += vI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i)
+          call pt2_add ( pt2_data_S(p),  1.d0,  pt2_data_teeth )
-          x3 += nI(pt2_stoch_istate, i) * pt2_W_T / pt2_w(i)
+          call pt2_add2( pt2_data_S2(p), 1.d0,  pt2_data_teeth )
-          Se(p) += x
-          Sv(p) += x2
-          Sn(p) += x3
-          Se2(p) += x*x
-          Sv2(p) += x2*x2
-          Sn2(p) += x3*3
         end do
-        avg  = E0 + Se(t) / dble(c)
+        call pt2_dealloc(pt2_data_teeth)
-        avg2 = v0 + Sv(t) / dble(c)
+        avg  = E0 + pt2_data_S(t) % pt2(pt2_stoch_istate) / dble(c)
-        avg3 = n0 + Sn(t) / dble(c)
+        avg2 = v0 + pt2_data_S(t) % variance(pt2_stoch_istate) / dble(c)
+        avg3 = n0 + pt2_data_S(t) % norm2(pt2_stoch_istate) / dble(c)
         if ((avg /= 0.d0) .or. (n == N_det_generators) ) then
           do_exit = .true.
         endif
@@ -500,18 +494,19 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
         call wall_time(time)
         ! 1/(N-1.5) : see  Brugger, The American Statistician (23) 4 p. 32 (1969)
         if(c > 2) then
-          eqt = dabs((Se2(t) / c) - (Se(t)/c)**2) ! dabs for numerical stability
+          eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
           eqt = sqrt(eqt / (dble(c) - 1.5d0))
           pt2_data % pt2_err(pt2_stoch_istate) = eqt
 
-          eqt = dabs((Sv2(t) / c) - (Sv(t)/c)**2) ! dabs for numerical stability
+          eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
           eqt = sqrt(eqt / (dble(c) - 1.5d0))
           pt2_data % variance_err(pt2_stoch_istate) = eqt
 
-          eqt = dabs((Sn2(t) / c) - (Sn(t)/c)**2) ! dabs for numerical stability
+          eqt = dabs((pt2_data_S2(t) % norm2(pt2_stoch_istate) / c) - (pt2_data_S(t) % norm2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
           eqt = sqrt(eqt / (dble(c) - 1.5d0))
           pt2_data % norm2_err(pt2_stoch_istate) = eqt
 
+
           if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
             time1 = time
             print '(G10.3, 2X, F16.10, 2X, G10.3, 2X, F14.10, 2X, F14.10, 2X, F10.4, A10)', c, avg+E, eqt, avg2, avg3, time-time0, ''
@@ -544,16 +539,22 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
           stop 'PT2: Unable to delete tasks (send)'
       endif
       do i=1,n_tasks
-        if(index(i).gt.size(eI,2).or.index(i).lt.1)then
+        if(index(i).gt.size(pt2_data_I,1).or.index(i).lt.1)then
          print*,'PB !!!'
          print*,'If you see this, send a bug report with the following content'
          print*,irp_here
-         print*,'i,index(i),size(ei,2) = ',i,index(i),size(ei,2)
+         print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1)
          stop -1
         endif
-        eI(1:N_states, index(i)) += pt2_data_task(i) % pt2(1:N_states)
+!        print *, pt2_data_I(index(i))%pt2
-        vI(1:N_states, index(i)) += pt2_data_task(i) % variance(1:N_states)
+!        print *, pt2_data_I(index(i))%variance
-        nI(1:N_states, index(i)) += pt2_data_task(i) % norm2(1:N_states)
+!        print *, pt2_data_I(index(i))%norm2
+!        print *, ''
+!        print *, pt2_data_task(i)%pt2
+!        print *, pt2_data_task(i)%variance
+!        print *, pt2_data_task(i)%norm2
+!        print *, ''
+        call pt2_add(pt2_data_I(index(i)),1.d0,pt2_data_task(i))
         f(index(i)) -= 1
       end do
       do i=1, b2%cur
@@ -566,6 +567,16 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, b, N_)
       endif
     end if
   end do
+  do i=1,N_det_generators
+    call pt2_dealloc(pt2_data_I(i))
+  enddo
+  do i=1,pt2_N_teeth+1
+    call pt2_dealloc(pt2_data_S(i))
+    call pt2_dealloc(pt2_data_S2(i))
+  enddo
+  do i=1,pt2_n_tasks_max
+    call pt2_dealloc(pt2_data_task(i))
+  enddo
 !print *,  'deleting b2'
   call delete_selection_buffer(b2)
 !print *,  'sorting b'

src/cipsi/pt2_type.irp.f

@@ -50,25 +50,51 @@ subroutine pt2_dealloc(pt2_data)
             )
 end subroutine
 
-subroutine pt2_add(p1, p2)
+subroutine pt2_add(p1, w, p2)
   implicit none
   use selection_types
   BEGIN_DOC
-! p1 += p2
+! p1 += w * p2
   END_DOC
   type(pt2_type), intent(inout) :: p1
+  double precision, intent(in)  :: w
   type(pt2_type), intent(in)    :: p2
 
-  p1 % pt2(:)            += p2 % pt2(:)
+  p1 % pt2(:)            = p1 % pt2(:)           + w * p2 % pt2(:)
-  p1 % pt2_err(:)        += p2 % pt2_err(:)
+  p1 % pt2_err(:)        = p1 % pt2_err(:)       + w * p2 % pt2_err(:)
-  p1 % rpt2(:)           += p2 % rpt2(:)
+  p1 % rpt2(:)           = p1 % rpt2(:)          + w * p2 % rpt2(:)
-  p1 % rpt2_err(:)       += p2 % rpt2_err(:)
+  p1 % rpt2_err(:)       = p1 % rpt2_err(:)      + w * p2 % rpt2_err(:)
-  p1 % variance(:)       += p2 % variance(:)
+  p1 % variance(:)       = p1 % variance(:)      + w * p2 % variance(:)
-  p1 % variance_err(:)   += p2 % variance_err(:)
+  p1 % variance_err(:)   = p1 % variance_err(:)  + w * p2 % variance_err(:)
-  p1 % norm2(:)          += p2 % norm2(:)
+  p1 % norm2(:)          = p1 % norm2(:)         + w * p2 % norm2(:)
-  p1 % norm2_err(:)      += p2 % norm2_err(:)
+  p1 % norm2_err(:)      = p1 % norm2_err(:)     + w * p2 % norm2_err(:)
-  p1 % overlap(:,:)      += p2 % overlap(:,:)
+  p1 % overlap(:,:)      = p1 % overlap(:,:)     + w * p2 % overlap(:,:)
-  p1 % overlap_err(:,:)  += p2 % overlap_err(:,:)
+  p1 % overlap_err(:,:)  = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:)
+
+end subroutine
+
+
+subroutine pt2_add2(p1, w, p2)
+  implicit none
+  use selection_types
+  BEGIN_DOC
+! p1 += w * p2**2
+  END_DOC
+  type(pt2_type), intent(inout) :: p1
+  double precision, intent(in)  :: w
+  type(pt2_type), intent(in)    :: p2
+
+  p1 % pt2(:)           = p1 % pt2(:)           + w * p2 % pt2(:)           * p2 % pt2(:)
+  p1 % pt2_err(:)       = p1 % pt2_err(:)       + w * p2 % pt2_err(:)       * p2 % pt2_err(:)
+  p1 % rpt2(:)          = p1 % rpt2(:)          + w * p2 % rpt2(:)          * p2 % rpt2(:)
+  p1 % rpt2_err(:)      = p1 % rpt2_err(:)      + w * p2 % rpt2_err(:)      * p2 % rpt2_err(:)
+  p1 % variance(:)      = p1 % variance(:)      + w * p2 % variance(:)      * p2 % variance(:)
+  p1 % variance_err(:)  = p1 % variance_err(:)  + w * p2 % variance_err(:)  * p2 % variance_err(:)
+  p1 % norm2(:)         = p1 % norm2(:)         + w * p2 % norm2(:)         * p2 % norm2(:)
+  p1 % norm2_err(:)     = p1 % norm2_err(:)     + w * p2 % norm2_err(:)     * p2 % norm2_err(:)
+  p1 % overlap(:,:)     = p1 % overlap(:,:)     + w * p2 % overlap(:,:)     * p2 % overlap(:,:)
+  p1 % overlap_err(:,:) = p1 % overlap_err(:,:) + w * p2 % overlap_err(:,:) * p2 % overlap_err(:,:)
+
 end subroutine
 
 

src/cipsi/run_pt2_slave.irp.f

@@ -462,7 +462,7 @@ subroutine pull_pt2_results(zmq_socket_pull, index, pt2_data, task_id, n_tasks,
  use f77_zmq
   implicit none
   integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
-  type(pt2_type), intent(inout)  :: pt2_data(n_tasks)
+  type(pt2_type), intent(inout)  :: pt2_data(*)
   type(selection_buffer), intent(inout) :: b
   integer, intent(out) :: index(*)
   integer, intent(out) :: n_tasks, task_id(*)