Fixed bug introduced by 0bf0513fb1

VERSION

@@ -0,0 +1 @@
+2.2.1

ocaml/qp_create_ezfio.ml

@@ -604,12 +604,33 @@ let run ?o b au c d m p cart xyz_file =
           ~rank:1 ~dim:[| shell_num |] ~data:ang_mom ) ;
         Ezfio.set_basis_shell_prim_index (Ezfio.ezfio_array_of_list
           ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_idx) ;
-        Ezfio.set_basis_shell_nucl (Ezfio.ezfio_array_of_list
+        Ezfio.set_basis_basis_nucleus_index (Ezfio.ezfio_array_of_list
-          ~rank:1 ~dim:[| shell_num |]
+          ~rank:1 ~dim:[| nucl_num |]
-          ~data:(list_map (fun (_,n) -> Nucl_number.to_int n) basis) ) ;
+          ~data:(
-        Ezfio.set_basis_shell_prim_coef (Ezfio.ezfio_array_of_list
+          list_map (fun (_,n) -> Nucl_number.to_int n) basis
+          |> List.fold_left (fun accu i -> 
+            match accu with 
+            | [] -> []
+            | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((h+1,i)::(h+1,j)::rest)
+          ) [(0,0)]
+          |> List.rev
+          |> List.map fst
+          )) ;
+        Ezfio.set_basis_nucleus_shell_num(Ezfio.ezfio_array_of_list
+          ~rank:1 ~dim:[| nucl_num |]
+          ~data:(
+          list_map (fun (_,n) -> Nucl_number.to_int n) basis
+          |> List.fold_left (fun accu i -> 
+            match accu with 
+            | [] -> [(1,i)]
+            | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((1,i)::(h,j)::rest)
+          ) []
+          |> List.rev
+          |> List.map fst
+          )) ;
+        Ezfio.set_basis_prim_coef (Ezfio.ezfio_array_of_list
           ~rank:1 ~dim:[| prim_num |] ~data:coef) ;
-        Ezfio.set_basis_shell_prim_expo (Ezfio.ezfio_array_of_list
+        Ezfio.set_basis_prim_expo (Ezfio.ezfio_array_of_list
           ~rank:1 ~dim:[| prim_num |] ~data:expo) ;
 
 

src/ao_basis/aos.irp.f

@@ -60,6 +60,10 @@ END_PROVIDER
       enddo
     endif
 
+    powA(1) = ao_power(i,1)
+    powA(2) = ao_power(i,2)
+    powA(3) = ao_power(i,3)
+
     ! Normalization of the contracted basis functions
     if (ao_normalized) then
       norm = 0.d0

src/basis/EZFIO.cfg

@@ -13,17 +13,18 @@ type: integer
 doc: Number of shells
 interface: ezfio, provider
 
+[nucleus_shell_num]
+type: integer
+doc: Number of shells per nucleus
+size: (nuclei.nucl_num)
+interface: ezfio, provider
+
 [shell_normalization_factor]
 type: double precision
-doc: Number of primitives per |AO|
+doc: Normalization factor applied to the whole shell, ex $1/\sqrt{ <d_{z^2}|d_{z^2}>}$
 size: (basis.shell_num)
 interface: ezfio
 
-[prim_num]
-type: integer
-doc: Total number of primitives
-interface: ezfio, provider
-
 [shell_ang_mom]
 type: integer
 doc: Angular momentum of each shell
@@ -42,19 +43,30 @@ doc: Max number of primitives in a shell
 size: (basis.shell_num)
 interface: ezfio, provider
 
-[shell_nucl]
+[basis_nucleus_index]
 type: integer
 doc: Index of the nucleus on which the shell is centered
-size: (basis.shell_num)
+size: (nuclei.nucl_num)
 interface: ezfio, provider
 
-[shell_prim_coef]
+[prim_normalization_factor]
+type: double precision
+doc: Normalization factor applied to each primitive
+size: (basis.prim_num)
+interface: ezfio
+
+[prim_num]
+type: integer
+doc: Total number of primitives
+interface: ezfio, provider
+
+[prim_coef]
 type: double precision
 doc: Primitive coefficients
 size: (basis.prim_num)
 interface: ezfio, provider
 
-[shell_prim_expo]
+[prim_expo]
 type: double precision
 doc: Exponents in the shell
 size: (basis.prim_num)

src/basis/basis.irp.f

@@ -15,7 +15,6 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
       call ezfio_get_basis_shell_normalization_factor(shell_normalization_factor)
     else
 
-
       double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
       integer                        :: l, powA(3), nz
       integer                        :: i,j,k
@@ -25,23 +24,22 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
       C_A(3) = 0.d0
 
       do i=1,shell_num
+
         powA(1) = shell_ang_mom(i)
         powA(2) = 0
         powA(3) = 0
 
-        ! Normalization of the contracted basis functions
         norm = 0.d0
-        do j=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1
+        do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
-          do k=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1
+          do j=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
-            call overlap_gaussian_xyz(C_A,C_A,shell_prim_expo(j),shell_prim_expo(k),&
+            call overlap_gaussian_xyz(C_A,C_A,prim_expo(j),prim_expo(k), &
-                powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+              powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
-            norm = norm+c*shell_prim_coef(j)*shell_prim_coef(k)
+            norm = norm+c*prim_coef(j)*prim_coef(k) * prim_normalization_factor(j) * prim_normalization_factor(k)
           enddo
         enddo
-        shell_normalization_factor(i) = dsqrt(norm)
+        shell_normalization_factor(i) = 1.d0/dsqrt(norm)
       enddo
 
-
     endif
   endif
   IRP_IF MPI_DEBUG
@@ -60,3 +58,62 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ]
   call write_time(6)
 
 END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Number of primitives per |AO|
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  if (mpi_master) then
+    if (size(prim_normalization_factor) == 0) return
+
+    call ezfio_has_basis_prim_normalization_factor(has)
+    if (has) then
+      write(6,'(A)') '.. >>>>> [ IO READ: prim_normalization_factor ] <<<<< ..'
+      call ezfio_get_basis_prim_normalization_factor(prim_normalization_factor)
+    else
+
+      double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
+      integer                        :: l, powA(3), nz
+      integer                        :: i,j,k
+      nz=100
+      C_A(1) = 0.d0
+      C_A(2) = 0.d0
+      C_A(3) = 0.d0
+
+      do i=1,shell_num
+
+        powA(1) = shell_ang_mom(i)
+        powA(2) = 0
+        powA(3) = 0
+
+        do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1
+            call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), &
+              powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+            prim_normalization_factor(k) = 1.d0/dsqrt(norm)
+        enddo
+      enddo
+
+
+    endif
+  endif
+  IRP_IF MPI_DEBUG
+  print *,  irp_here, mpi_rank
+  call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  IRP_IF MPI
+  include 'mpif.h'
+  integer                        :: ierr
+  call MPI_BCAST( prim_normalization_factor, (prim_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+  if (ierr /= MPI_SUCCESS) then
+    stop 'Unable to read prim_normalization_factor with MPI'
+  endif
+  IRP_ENDIF
+
+  call write_time(6)
+
+END_PROVIDER

src/zmq/utils.irp.f

@@ -127,7 +127,9 @@ function zmq_port(ishift)
   END_DOC
   integer, intent(in)            :: ishift
   character*(8)                  :: zmq_port
+  !$OMP CRITICAL(write)
   write(zmq_port,'(I8)') zmq_port_start+ishift
+  !$OMP END CRITICAL(write)
   zmq_port = adjustl(trim(zmq_port))
 end
 
@@ -518,7 +520,9 @@ subroutine new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
   zmq_socket_pull      = new_zmq_pull_socket ()
+  !$OMP CRITICAL(write)
   write(name,'(A,I8.8)') trim(name_in)//'.', icount
+  !$OMP END CRITICAL(write)
   sze = len(trim(name))
   zmq_state = trim(name)
   call lowercase(name,sze)
@@ -582,7 +586,9 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
   integer, save                  :: icount=0
 
   icount = icount+1
+  !$OMP CRITICAL(write)
   write(name,'(A,I8.8)') trim(name_in)//'.', icount
+  !$OMP END CRITICAL(write)
   sze = len(trim(name))
   call lowercase(name,sze)
   if (name /= zmq_state) then
@@ -704,7 +710,9 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i
 
   disconnect_from_taskserver_state = -1
 
+  !$OMP CRITICAL(write)
   write(message,*) 'disconnect '//trim(state), worker_id
+  !$OMP END CRITICAL(write)
 
   sze = min(510,len(trim(message)))
   rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@@ -781,7 +789,9 @@ integer function zmq_abort(zmq_to_qp_run_socket)
   character*(512)                :: message
   zmq_abort = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'abort '
+  !$OMP END CRITICAL(write)
 
 
   sze = len(trim(message))
@@ -823,7 +833,9 @@ integer function task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_i
 
   task_done_to_taskserver = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'task_done '//trim(zmq_state), worker_id, task_id
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@@ -856,9 +868,11 @@ integer function tasks_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_
 
   tasks_done_to_taskserver = 0
 
+  !$OMP CRITICAL(write)
   allocate(character(LEN=64+n_tasks*12) :: message)
   write(fmt,*) '(A,X,A,I10,X,', n_tasks, '(I11,1X))'
   write(message,*) 'task_done '//trim(zmq_state), worker_id, (task_id(k), k=1,n_tasks)
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@@ -900,7 +914,9 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id
 
   get_task_from_taskserver = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'get_task '//trim(zmq_state), worker_id
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
@@ -961,7 +977,9 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i
 
   get_tasks_from_taskserver = 0
 
+  !$OMP CRITICAL(write)
   write(message,'(A,A,X,I10,I10)') 'get_tasks ', trim(zmq_state), worker_id, n_tasks
+  !$OMP END CRITICAL(write)
 
   sze = len(trim(message))
   rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
@@ -1061,7 +1079,9 @@ integer function zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,mo
 
   zmq_delete_task = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'del_task ', zmq_state, task_id
+  !$OMP END CRITICAL(write)
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
   if (rc /= len(trim(message))) then
     zmq_delete_task = -1
@@ -1101,7 +1121,9 @@ integer function zmq_delete_task_async_send(zmq_to_qp_run_socket,task_id,sending
   endif
   zmq_delete_task_async_send = 0
 
+  !$OMP CRITICAL(write)
   write(message,*) 'del_task ', zmq_state, task_id
+  !$OMP END CRITICAL(write)
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
   if (rc /= len(trim(message))) then
     zmq_delete_task_async_send = -1
@@ -1159,8 +1181,10 @@ integer function zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,n
 
   allocate(character(LEN=64+n_tasks*12) :: message)
 
+  !$OMP CRITICAL(write)
   write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))'
   write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks)
+  !$OMP END CRITICAL(write)
 
 
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
@@ -1206,8 +1230,10 @@ integer function zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_task
 
   allocate(character(LEN=64+n_tasks*12) :: message)
 
+  !$OMP CRITICAL(write)
   write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))'
   write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks)
+  !$OMP END CRITICAL(write)
 
 
   rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)