diff --git a/plugins/Symmetry/NEEDED_CHILDREN_MODULES b/plugins/Symmetry/NEEDED_CHILDREN_MODULES
new file mode 100644
index 00000000..7ea84ba4
--- /dev/null
+++ b/plugins/Symmetry/NEEDED_CHILDREN_MODULES
@@ -0,0 +1 @@
+Bitmask Nuclei Determinants
diff --git a/plugins/Symmetry/README.rst b/plugins/Symmetry/README.rst
new file mode 100644
index 00000000..ba643f88
--- /dev/null
+++ b/plugins/Symmetry/README.rst
@@ -0,0 +1,12 @@
+========
+Symmetry
+========
+
+Needed Modules
+==============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
+Documentation
+=============
+.. Do not edit this section It was auto-generated
+.. by the `update_README.py` script.
diff --git a/plugins/Symmetry/Symmetry.main.irp.f b/plugins/Symmetry/Symmetry.main.irp.f
new file mode 100644
index 00000000..d8d2f350
--- /dev/null
+++ b/plugins/Symmetry/Symmetry.main.irp.f
@@ -0,0 +1,30 @@
+program Symmetry
+  implicit none
+  BEGIN_DOC
+! TODO
+  END_DOC
+  integer :: n_sym_points, i, j
+  double precision, allocatable :: tmp(:,:), sym_points(:,:)
+  character*8 :: sym
+  n_sym_points = 3
+
+  print *,  'Molecule is linear:', molecule_is_linear
+  print *,  'Has center of inversion:', molecule_has_center_of_inversion
+  print *,  'Symmetry rotation axis:', sym_rotation_axis(:)
+  call find_symmetry(sym)
+  print *,  'Group:'//sym
+  return
+  allocate(tmp(n_sym_points,n_det), sym_points(3,n_sym_points))
+
+  call generate_sym_coord(n_sym_points,sym_points)
+  call compute_sym_det_values(sym_points,n_sym_points,tmp)
+  do i=1,mo_tot_num
+    print *,  i, real(tmp(1:3,i))
+  enddo
+  sym_points(1,:) = -sym_points(1,:)
+  call compute_sym_det_values(sym_points,n_sym_points,tmp)
+  do i=1,n_det
+    print *,  i, real(tmp(1:3,i))
+  enddo
+
+end
diff --git a/plugins/Symmetry/aos.irp.f b/plugins/Symmetry/aos.irp.f
new file mode 100644
index 00000000..97b7aa62
--- /dev/null
+++ b/plugins/Symmetry/aos.irp.f
@@ -0,0 +1,113 @@
+BEGIN_PROVIDER [ double precision, sym_box, (3,2) ]
+  implicit none
+  BEGIN_DOC
+  ! Opposite points of the box containing the molecule
+  END_DOC
+  integer                        :: i,xyz
+  sym_box(:,:) = 0.d0
+  do xyz=1,3
+    do i=1,nucl_num
+      sym_box(xyz,1) = min(sym_box(xyz,1), nucl_coord(i,xyz))
+      sym_box(xyz,2) = max(sym_box(xyz,2), nucl_coord(i,xyz))
+    enddo
+  enddo
+  sym_box(:,1) = sym_box(:,1) - 2.d0
+  sym_box(:,2) = sym_box(:,2) + 2.d0
+END_PROVIDER
+
+subroutine generate_sym_coord(n_sym_points,result)
+  implicit none
+  integer, intent(in)            :: n_sym_points
+  double precision, intent(out)  :: result(3,n_sym_points)
+  BEGIN_DOC
+  ! xyz coordinates of points to check the symmetry, drawn uniformly in the molecular box.
+  END_DOC
+  integer                        :: i, xyz
+  
+  call random_number(result)
+  do xyz=1,3
+    result(xyz,:) = sym_box(xyz,1) + result(xyz,:) * (sym_box(xyz,2)-sym_box(xyz,1))
+  enddo
+  
+end
+
+
+subroutine compute_sym_ao_values(sym_points, n_sym_points, result)
+  implicit none
+  BEGIN_DOC
+  ! Values of the AO symmetry functions
+  END_DOC
+  integer, intent(in)            :: n_sym_points
+  double precision, intent(in)   :: sym_points(3,n_sym_points)
+  double precision, intent(out)  :: result(n_sym_points, ao_num)
+  integer                        :: i, j
+  double precision               :: x, y, z
+  result (:,:) = 0.d0
+  do j=1,ao_num
+    do i=1,n_sym_points
+      x = sym_points(1,i) - nucl_coord_transp(1,ao_nucl(j))
+      y = sym_points(2,i) - nucl_coord_transp(2,ao_nucl(j))
+      z = sym_points(3,i) - nucl_coord_transp(3,ao_nucl(j))
+      x = x**ao_power(j,1)
+      y = y**ao_power(j,2)
+      z = z**ao_power(j,3)
+      result(i,j) = x*y*z*exp(-(x*x+y*y+z*z))
+    enddo
+  enddo
+  
+end
+
+subroutine compute_sym_mo_values(sym_points, n_sym_points, result)
+  implicit none
+  BEGIN_DOC
+  ! Values of the MO symmetry functions
+  END_DOC
+  integer, intent(in)            :: n_sym_points
+  double precision, intent(in)   :: sym_points(3,n_sym_points)
+  double precision, intent(out)  :: result(n_sym_points, mo_tot_num)
+
+  double precision, allocatable :: tmp(:,:)
+  allocate(tmp(n_sym_points,ao_num))
+  call compute_sym_ao_values(sym_points,n_sym_points,tmp)
+  call dgemm('N','N',n_sym_points,mo_tot_num,ao_num, &
+    1.d0, tmp,size(tmp,1), mo_coef, size(mo_coef,1), &
+    0.d0, result,size(result,1))
+  deallocate(tmp)
+end
+
+
+subroutine compute_sym_det_values(sym_points, n_sym_points, result)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Values of the determinant symmetry functions
+  END_DOC
+  integer, intent(in)            :: n_sym_points
+  double precision, intent(in)   :: sym_points(3,n_sym_points)
+  double precision, intent(out)  :: result(n_sym_points, N_det)
+
+  integer                        :: list(N_int*bit_kind_size,2)
+  integer                        :: n_elements(2)
+
+  integer :: i, j, imo
+
+  double precision, allocatable :: tmp(:,:)
+
+  allocate(tmp(n_sym_points,mo_tot_num))
+  call compute_sym_mo_values(sym_points, n_sym_points, tmp)
+
+  result = 1.d0
+  do i=1,N_det
+    call bitstring_to_list_ab(psi_det(1,1,i), list, n_elements, N_int)
+    do j=1,n_elements(1)
+      imo = list(j,1)
+      result(:,i) *= tmp(:,imo)
+    enddo
+    do j=1,n_elements(2)
+      imo = list(j,2)
+      result(:,i) *= tmp(:,imo)
+    enddo
+  enddo
+
+  deallocate(tmp)
+end
diff --git a/plugins/Symmetry/find_sym.irp.f b/plugins/Symmetry/find_sym.irp.f
new file mode 100644
index 00000000..4a100805
--- /dev/null
+++ b/plugins/Symmetry/find_sym.irp.f
@@ -0,0 +1,318 @@
+BEGIN_PROVIDER [ logical, molecule_is_linear ]
+ implicit none
+ BEGIN_DOC
+ ! True if the molecule is linear
+ END_DOC
+ molecule_is_linear = (minval(inertia_tensor_eigenvalues) < 1.d-5)
+END_PROVIDER
+
+BEGIN_PROVIDER [ logical, molecule_has_center_of_inversion ]
+  implicit none
+  BEGIN_DOC
+  ! If true, there is a center of inversion in the WF
+  END_DOC
+  molecule_has_center_of_inversion = .True.
+  integer                        :: i,j,k
+  double precision               :: point(3)
+  logical                        :: found
+  double precision, external     :: u_dot_u
+  do i=1,nucl_num
+    found = .False.
+    do j=1,nucl_num
+      if (nucl_charge(i) /= nucl_charge(j)) cycle
+      point(:) = nucl_coord_transp(:,i) + nucl_coord_transp(:,j)
+      if (u_dot_u(point,3) < 1.d-5) then
+        found = .True.
+        exit
+      endif
+    enddo
+    if (.not.found) then
+      molecule_has_center_of_inversion = .False.
+      exit
+    endif
+  enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, sym_rotation_axis, (3) ]
+  implicit none
+  BEGIN_DOC
+  ! Order of the rotation axis
+  END_DOC
+  include 'constants.include.F'
+  
+  integer                        :: i,j,k
+  double precision               :: point(3), point2(3)
+  logical                        :: found
+  double precision, external     :: u_dot_u
+  integer                        :: iorder, iaxis
+  double precision               :: theta, sin_t, cos_t
+  do iaxis=1,3
+    do iorder=12,2,-1
+      sym_rotation_axis(iaxis) = iorder
+      theta = 2.d0*pi/dble(iorder)
+      sin_t = dsin(theta)
+      cos_t = dcos(theta)
+      do i=1,nucl_num
+        found = .False.
+        if (iaxis==1) then
+          point(1) = nucl_coord_transp(1,i)
+          point(2) = nucl_coord_transp(2,i)*cos_t - nucl_coord_transp(3,i)*sin_t
+          point(3) = nucl_coord_transp(2,i)*sin_t + nucl_coord_transp(3,i)*cos_t
+        else if (iaxis==2) then
+          point(1) = nucl_coord_transp(1,i)*cos_t - nucl_coord_transp(3,i)*sin_t
+          point(2) = nucl_coord_transp(2,i)
+          point(3) = nucl_coord_transp(1,i)*sin_t + nucl_coord_transp(3,i)*cos_t
+        endif
+        if (iaxis==3) then
+          point(1) = nucl_coord_transp(1,i)*cos_t - nucl_coord_transp(2,i)*sin_t
+          point(2) = nucl_coord_transp(1,i)*sin_t + nucl_coord_transp(2,i)*cos_t
+          point(3) = nucl_coord_transp(3,i)
+        endif
+        do j=1,nucl_num
+          if (nucl_charge(i) /= nucl_charge(j)) cycle
+          point2(:) = nucl_coord_transp(:,j) - point(:)
+          if (u_dot_u(point2,3) < 1.d-5) then
+            found = .True.
+            exit
+          endif
+        enddo
+        if (.not.found) then
+          sym_rotation_axis(iaxis) = 1
+          exit
+        endif
+      enddo
+      if (sym_rotation_axis(iaxis) /= 1) then
+        exit
+      endif
+    enddo
+  enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, molecule_principal_axis ]
+&BEGIN_PROVIDER [ integer, molecule_secondary_axis ]
+&BEGIN_PROVIDER [ integer, molecule_ternary_axis ]
+  implicit none
+  BEGIN_DOC
+! Which axis is the Z axis
+  END_DOC
+  molecule_principal_axis = maxloc(sym_rotation_axis,1)
+  if (molecule_principal_axis == 1) then
+    if (sym_rotation_axis(2) >= sym_rotation_axis(3)) then 
+      molecule_secondary_axis = 2
+      molecule_ternary_axis = 3
+    else
+      molecule_secondary_axis = 3
+      molecule_ternary_axis = 2
+    endif
+  else if (molecule_principal_axis == 2) then
+    if (sym_rotation_axis(1) >= sym_rotation_axis(3)) then 
+      molecule_secondary_axis = 1
+      molecule_ternary_axis = 3
+    else
+      molecule_secondary_axis = 3
+      molecule_ternary_axis = 1
+    endif
+  else if (molecule_principal_axis == 3) then
+    if (sym_rotation_axis(1) >= sym_rotation_axis(2)) then 
+      molecule_secondary_axis = 1
+      molecule_ternary_axis = 2
+    else
+      molecule_secondary_axis = 2
+      molecule_ternary_axis = 1
+    endif
+  endif
+END_PROVIDER
+
+BEGIN_PROVIDER [ logical, molecule_has_secondary_c2_rotation ]
+  implicit none
+  BEGIN_DOC
+! If true, the molecule has a C2 rotation perpendicular to the main axis
+  END_DOC
+  if (molecule_principal_axis == 1) then
+    molecule_has_secondary_c2_rotation = (sym_rotation_axis(2)==2) .or. (sym_rotation_axis(3)==2)
+  else if (molecule_principal_axis == 2) then
+    molecule_has_secondary_c2_rotation = (sym_rotation_axis(1)==2) .or. (sym_rotation_axis(3)==2)
+  else if (molecule_principal_axis == 3) then
+    molecule_has_secondary_c2_rotation = (sym_rotation_axis(1)==2) .or. (sym_rotation_axis(2)==2)
+  endif
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ logical, molecule_has_improper_rotation ]
+  implicit none
+  BEGIN_DOC
+  ! Order of the rotation axis
+  END_DOC
+  include 'constants.include.F'
+  
+  integer                        :: i,j,k
+  double precision               :: point(3), point2(3)
+  logical                        :: found
+  double precision, external     :: u_dot_u
+  integer                        :: iorder, iaxis
+  double precision               :: theta, sin_t, cos_t
+  iaxis=molecule_principal_axis
+  iorder = 2*sym_rotation_axis(iaxis)
+  molecule_has_improper_rotation = .True.
+  theta = 2.d0*pi/dble(iorder)
+  sin_t = dsin(theta)
+  cos_t = dcos(theta)
+  do i=1,nucl_num
+    found = .False.
+    if (iaxis==1) then
+      point(1) = -nucl_coord_transp(1,i)
+      point(2) = nucl_coord_transp(2,i)*cos_t - nucl_coord_transp(3,i)*sin_t
+      point(3) = nucl_coord_transp(2,i)*sin_t + nucl_coord_transp(3,i)*cos_t
+    else if (iaxis==2) then
+      point(1) = nucl_coord_transp(1,i)*cos_t - nucl_coord_transp(3,i)*sin_t
+      point(2) = -nucl_coord_transp(2,i)
+      point(3) = nucl_coord_transp(1,i)*sin_t + nucl_coord_transp(3,i)*cos_t
+    endif
+    if (iaxis==3) then
+      point(1) = nucl_coord_transp(1,i)*cos_t - nucl_coord_transp(2,i)*sin_t
+      point(2) = nucl_coord_transp(1,i)*sin_t + nucl_coord_transp(2,i)*cos_t
+      point(3) = -nucl_coord_transp(3,i)
+    endif
+    do j=1,nucl_num
+      if (nucl_charge(i) /= nucl_charge(j)) cycle
+      point2(:) = nucl_coord_transp(:,j) - point(:)
+      if (u_dot_u(point2,3) < 1.d-5) then
+        found = .True.
+        exit
+      endif
+    enddo
+    if (.not.found) then
+      molecule_has_improper_rotation = .False.
+      exit
+    endif
+  enddo
+  
+END_PROVIDER
+
+BEGIN_PROVIDER [ logical, molecule_has_center_of_inversion ]
+  implicit none
+  BEGIN_DOC
+  ! If true, there is a center of inversion in the WF
+  END_DOC
+  molecule_has_center_of_inversion = .True.
+  integer                        :: i,j,k
+  double precision               :: point(3)
+  logical                        :: found
+  double precision, external     :: u_dot_u
+  do i=1,nucl_num
+    found = .False.
+    do j=1,nucl_num
+      if (nucl_charge(i) /= nucl_charge(j)) cycle
+      point(:) = nucl_coord_transp(:,i) + nucl_coord_transp(:,j)
+      if (u_dot_u(point,3) < 1.d-5) then
+        found = .True.
+        exit
+      endif
+    enddo
+    if (.not.found) then
+      molecule_has_center_of_inversion = .False.
+      exit
+    endif
+  enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ logical, molecule_has_sigma_plane, (3) ]
+  implicit none
+  BEGIN_DOC
+  ! If true, there is a symmetry plane perpendicular to the main axis
+  END_DOC
+  integer                        :: i,j,k
+  double precision               :: point(3), point2(3)
+  logical                        :: found
+  double precision, external     :: u_dot_u
+  integer                        :: iaxis
+  do iaxis=1,3
+    molecule_has_sigma_plane(iaxis) = .True.
+    do i=1,nucl_num
+      found = .False.
+      point(:) = nucl_coord_transp(:,i)
+      point(iaxis) = -point(iaxis)
+      do j=1,nucl_num
+        if (nucl_charge(i) /= nucl_charge(j)) cycle
+        point2(:) = nucl_coord_transp(:,j) - point(:)
+        if (u_dot_u(point2,3) < 1.d-5) then
+          found = .True.
+          exit
+        endif
+      enddo
+      if (.not.found) then
+        molecule_has_sigma_plane(iaxis) = .False.
+        exit
+      endif
+    enddo
+  enddo
+
+END_PROVIDER
+
+subroutine find_symmetry(result)
+  implicit none
+  BEGIN_DOC
+! Finds the point group of the molecule
+  END_DOC
+  character*16       :: result
+  character*2, save  :: i_to_a(24) = (/ '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12',&
+      '13', '14', '15', '16', '17', '18', '19', '20', '21', '22', '23', '24'  /)
+  result = 'C1'
+  if (molecule_is_linear) then
+    if (molecule_has_center_of_inversion) then
+      result = 'D*h'
+    else
+      result = 'C*v'
+    endif
+  else
+    if (maxval(sym_rotation_axis) == 1) then
+      if (molecule_has_sigma_plane(1).or.molecule_has_sigma_plane(2).or.&
+            molecule_has_sigma_plane(3) ) then
+        result = 'Cs'
+      else
+        if (molecule_has_center_of_inversion) then
+          result = 'Ci'
+        endif
+      endif
+    else
+      if (molecule_has_secondary_c2_rotation) then
+        if (molecule_has_sigma_plane(molecule_principal_axis)) then
+          result = 'D'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))//'h'
+        else
+          if (molecule_has_sigma_plane(molecule_secondary_axis)) then
+            result = 'D'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))//'d'
+          else
+            if (molecule_has_sigma_plane(molecule_ternary_axis)) then
+              result = 'D'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))//'d'
+            else
+              result = 'D'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))
+            endif
+          endif
+        endif
+      else
+        if (molecule_has_sigma_plane(molecule_principal_axis)) then
+          result = 'C'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))//'h'
+        else
+          if (molecule_has_sigma_plane(molecule_secondary_axis)) then
+            result = 'C'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))//'v'
+          else
+            if (molecule_has_sigma_plane(molecule_ternary_axis)) then
+              result = 'C'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))//'v'
+            else
+              if (molecule_has_improper_rotation) then
+                result = 'S'//trim(i_to_a(2*sym_rotation_axis(molecule_principal_axis)))
+              else
+                result = 'C'//trim(i_to_a(sym_rotation_axis(molecule_principal_axis)))
+              endif
+            endif
+          endif
+        endif
+      endif
+    endif
+  endif
+end
diff --git a/src/Nuclei/inertia.irp.f b/src/Nuclei/inertia.irp.f
new file mode 100644
index 00000000..c4517cde
--- /dev/null
+++ b/src/Nuclei/inertia.irp.f
@@ -0,0 +1,32 @@
+BEGIN_PROVIDER [ double precision, inertia_tensor, (3,3) ]
+  implicit none
+  BEGIN_DOC
+  ! Inertia tensor
+  END_DOC
+  integer                        :: i,j,k
+  inertia_tensor = 0.d0
+  do k=1,nucl_num
+    inertia_tensor(1,1) += mass(int(nucl_charge(k))) * ((nucl_coord_input(k,2)-center_of_mass(2))**2 + (nucl_coord_input(k,3)-center_of_mass(3))**2)
+    inertia_tensor(2,2) += mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))**2 + (nucl_coord_input(k,3)-center_of_mass(3))**2)
+    inertia_tensor(3,3) += mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))**2 + (nucl_coord_input(k,2)-center_of_mass(2))**2)
+    inertia_tensor(1,2) -= mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))    * (nucl_coord_input(k,2)-center_of_mass(2))   )
+    inertia_tensor(1,3) -= mass(int(nucl_charge(k))) * ((nucl_coord_input(k,1)-center_of_mass(1))    * (nucl_coord_input(k,3)-center_of_mass(3))   )
+    inertia_tensor(2,3) -= mass(int(nucl_charge(k))) * ((nucl_coord_input(k,2)-center_of_mass(2))    * (nucl_coord_input(k,3)-center_of_mass(3))   )
+  enddo
+  inertia_tensor(2,1) = inertia_tensor(1,2) 
+  inertia_tensor(3,1) = inertia_tensor(1,3) 
+  inertia_tensor(3,2) = inertia_tensor(2,3) 
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, inertia_tensor_eigenvectors, (3,3) ]
+&BEGIN_PROVIDER [ double precision, inertia_tensor_eigenvalues , (3) ]
+ implicit none
+ BEGIN_DOC
+ ! Eigenvectors/eigenvalues of the inertia_tensor. Used to find normal orientation.
+ END_DOC
+ integer :: k
+ call lapack_diagd(inertia_tensor_eigenvalues,inertia_tensor_eigenvectors,inertia_tensor,3,3)
+ print *, 'Rotational constants (GHZ):'
+ print *, (1805.65468542d0/inertia_tensor_eigenvalues(k), k=3,1,-1) 
+END_PROVIDER
+
diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f
index 307a9610..43609506 100644
--- a/src/Nuclei/nuclei.irp.f
+++ b/src/Nuclei/nuclei.irp.f
@@ -1,4 +1,4 @@
-BEGIN_PROVIDER [ double precision, nucl_coord,  (nucl_num,3) ]
+BEGIN_PROVIDER [ double precision, nucl_coord_input,  (nucl_num,3) ]
    implicit none
    
    BEGIN_DOC
@@ -8,7 +8,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord,  (nucl_num,3) ]
    
    if (mpi_master) then
      double precision, allocatable  :: buffer(:,:)
-     nucl_coord = 0.d0
+     nucl_coord_input = 0.d0
      allocate (buffer(nucl_num,3))
      buffer = 0.d0
      logical                        :: has
@@ -22,7 +22,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord,  (nucl_num,3) ]
      
      do i=1,3
        do j=1,nucl_num
-         nucl_coord(j,i) = buffer(j,i)
+         nucl_coord_input(j,i) = buffer(j,i)
        enddo
      enddo
      deallocate(buffer)
@@ -31,6 +31,65 @@ BEGIN_PROVIDER [ double precision, nucl_coord,  (nucl_num,3) ]
      character*(64), parameter      :: ft= '(A16, 4(1X,A12  ))'
      double precision, parameter    :: a0= 0.529177249d0
      
+     call write_time(output_Nuclei)
+     write(output_Nuclei,'(A)') ''
+     write(output_Nuclei,'(A)') 'Input Nuclear Coordinates (Angstroms)'
+     write(output_Nuclei,'(A)') '====================================='
+     write(output_Nuclei,'(A)') ''
+     write(output_Nuclei,ft)                                         &
+         '================','============','============','============','============'
+     write(output_Nuclei,*)                                          &
+         '     Atom          Charge          X            Y            Z '
+     write(output_Nuclei,ft)                                         &
+         '================','============','============','============','============'
+     do i=1,nucl_num
+       write(output_Nuclei,f) nucl_label(i), nucl_charge(i),         &
+           nucl_coord_input(i,1)*a0,                                       &
+           nucl_coord_input(i,2)*a0,                                       &
+           nucl_coord_input(i,3)*a0
+     enddo
+     write(output_Nuclei,ft)                                         &
+         '================','============','============','============','============'
+     write(output_Nuclei,'(A)') ''
+     
+   endif
+   
+   IRP_IF MPI
+     include 'mpif.h'
+     integer                        :: ierr
+     call MPI_BCAST( nucl_coord_input, 3*nucl_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+     if (ierr /= MPI_SUCCESS) then
+       stop 'Unable to read nucl_coord_input with MPI'
+     endif
+   IRP_ENDIF
+
+END_PROVIDER
+ 
+BEGIN_PROVIDER [ double precision, nucl_coord,  (nucl_num,3) ]
+   implicit none
+   
+   BEGIN_DOC
+   ! Nuclear coordinates in standard orientation
+   END_DOC
+   
+   if (mpi_master) then
+     integer :: i
+     do i=1,nucl_num
+       nucl_coord(i,1) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,1) + &
+                         (nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,1) + &
+                         (nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,1) 
+       nucl_coord(i,2) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,2) + &
+                         (nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,2) + &
+                         (nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,2) 
+       nucl_coord(i,3) = (nucl_coord_input(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,3) + &
+                         (nucl_coord_input(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,3) + &
+                         (nucl_coord_input(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,3) 
+     enddo
+     
+     character*(64), parameter      :: f = '(A16, 4(1X,F12.6))'
+     character*(64), parameter      :: ft= '(A16, 4(1X,A12  ))'
+     double precision, parameter    :: a0= 0.529177249d0
+     
      call write_time(output_Nuclei)
      write(output_Nuclei,'(A)') ''
      write(output_Nuclei,'(A)') 'Nuclear Coordinates (Angstroms)'
@@ -285,3 +344,87 @@ BEGIN_PROVIDER [ character*(128), element_name, (78)]
  element_name(78) = 'Pt'
 
 END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mass, (0:110) ]
+  implicit none
+  BEGIN_DOC
+  ! Atomic masses
+  END_DOC
+  
+  mass(  0   )  =  0.
+  mass(  1   )  =  1.0079
+  mass(  2   )  =  4.00260
+  mass(  3   )  =  6.941
+  mass(  4   )  =  9.01218
+  mass(  5   )  =  10.81
+  mass(  6   )  =  12.011
+  mass(  7   )  =  14.0067
+  mass(  8   )  =  15.9994
+  mass(  9   )  =  18.998403
+  mass(  10  )  =  20.179
+  mass(  11  )  =  22.98977
+  mass(  12  )  =  24.305
+  mass(  13  )  =  26.98154
+  mass(  14  )  =  28.0855
+  mass(  15  )  =  30.97376
+  mass(  16  )  =  32.06
+  mass(  17  )  =  35.453
+  mass(  18  )  =  39.948
+  mass(  19  )  =  39.0983
+  mass(  20  )  =  40.08
+  mass(  21  )  =  44.9559
+  mass(  22  )  =  47.90
+  mass(  23  )  =  50.9415
+  mass(  24  )  =  51.996
+  mass(  25  )  =  54.9380
+  mass(  26  )  =  55.9332
+  mass(  27  )  =  58.9332
+  mass(  28  )  =  58.70
+  mass(  29  )  =  63.546
+  mass(  30  )  =  65.38
+  mass(  31  )  =  69.72
+  mass(  32  )  =  72.59
+  mass(  33  )  =  74.9216
+  mass(  34  )  =  78.96
+  mass(  35  )  =  79.904
+  mass(  36  )  =  83.80
+  mass(  37  )  =  85.4678
+  mass(  38  )  =  87.62
+  mass(  39  )  =  88.90584
+  mass(  40  )  =  91.224
+  mass(  41  )  =  92.90637
+  mass(  42  )  =  95.95
+  mass(  43  )  =  98.
+  mass(  44  )  =  101.07
+  mass(  45  )  =  102.90550
+  mass(  46  )  =  106.42
+  mass(  47  )  =  107.8682
+  mass(  48  )  =  112.414
+  mass(  49  )  =  114.818
+  mass(  50  )  =  118.710
+  mass(  51  )  =  121.760
+  mass(  52  )  =  127.60
+  mass(  53  )  =  126.90447
+  mass(  54  )  =  131.293
+  mass(  78  )  =  195.084
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, center_of_mass, (3) ]
+  implicit none
+  BEGIN_DOC
+  ! Center of mass of the molecule
+  END_DOC
+  integer                        :: i,j
+  double precision               :: s
+  center_of_mass(:) = 0.d0
+  s = 0.d0
+  do i=1,nucl_num
+    do j=1,3
+      center_of_mass(j) += nucl_coord_input(i,j)* mass(int(nucl_charge(i)))
+    enddo
+    s += mass(int(nucl_charge(i)))
+  enddo
+  s = 1.d0/s
+  center_of_mass(:) = center_of_mass(:)*s
+END_PROVIDER
+