program print_mos
 implicit none
 character*(128) :: output
 integer :: i_unit_output,getUnitAndOpen
 provide ezfio_filename 

 integer :: i
 print*,trim(ezfio_filename)
 output=trim(ezfio_filename)//'.mol'
 print*,'output = ',trim(output)
 i_unit_output = getUnitAndOpen(output,'w')
 print*,'i_unit_output = ',i_unit_output
 call write_intro_gamess(i_unit_output)
 call write_geometry(i_unit_output)
 call write_Ao_basis(i_unit_output)
 call write_Mo_basis(i_unit_output)
 

 write(i_unit_output,*)''
 write(i_unit_output,*)''
 write(i_unit_output,*)'          ------------------------'

 close(i_unit_output)
end

subroutine write_intro_gamess(i_unit_output)
 implicit none
 integer, intent(in) :: i_unit_output
 integer :: i,j,k,l

 write(i_unit_output,*)'         *         GAMESS VERSION = 22 FEB 2006 (R5)          *' 
 write(i_unit_output,*)'         *             FROM IOWA STATE UNIVERSITY             *' 
 write(i_unit_output,*)'         * M.W.SCHMIDT, K.K.BALDRIDGE, J.A.BOATZ, S.T.ELBERT, *' 
 write(i_unit_output,*)'         *   M.S.GORDON, J.H.JENSEN, S.KOSEKI, N.MATSUNAGA,   *' 
 write(i_unit_output,*)'         *          K.A.NGUYEN, S.J.SU, T.L.WINDUS,           *' 
 write(i_unit_output,*)'         *       TOGETHER WITH M.DUPUIS, J.A.MONTGOMERY       *' 
 write(i_unit_output,*)'         *         J.COMPUT.CHEM.  14, 1347-1363(1993)        *' 
 write(i_unit_output,*)''

end




subroutine write_geometry(i_unit_output)
 implicit none
 integer, intent(in) :: i_unit_output
 integer :: i,j,k,l, getUnitAndOpen



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 write(i_unit_output,*)'ATOM      ATOMIC                      COORDINATES (BOHR)          '
 write(i_unit_output,*)'          CHARGE         X                   Y                   Z'
 do i = 1, nucl_num
! write(i_unit_output,'(A2 I3 X F3.1 X 3(F16.10))') trim(element_name(int(nucl_charge(i)))),i,(nucl_charge(i)), nucl_coord(i,1), nucl_coord(i,2), nucl_coord(i,3)
  write(i_unit_output,'(A2,I1, 9X F5.1 X 3(F16.10 ,4X))') trim(element_name(int(nucl_charge(i)))),i,(nucl_charge(i)), nucl_coord(i,1), nucl_coord(i,2), nucl_coord(i,3)
 enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end


subroutine write_Ao_basis(i_unit_output)
 implicit none
 integer, intent(in) :: i_unit_output
 integer :: i,j,k,l, getUnitAndOpen
 character*(128) :: character_shell
 integer :: i_shell,i_prim,i_ao
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 write(i_unit_output,*)''
 write(i_unit_output,*)''
 write(i_unit_output,*)'    ATOMIC BASIS SET'
 write(i_unit_output,*)'    ----------------'
 write(i_unit_output,*)'THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN UNNORMALIZED'
 write(i_unit_output,*)'THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED TO UNITY'
 write(i_unit_output,*)''
 write(i_unit_output,*)'SHELL TYPE  PRIMITIVE        EXPONENT          CONTRACTION COEFFICIENT(S)'
 write(i_unit_output,*)''
 write(i_unit_output,*)''

 i_shell = 0
 i_prim = 0
 do i = 1, Nucl_num
   write(i_unit_output,'(A2,I1)') trim(element_name(int(nucl_charge(i)))),i
   write(i_unit_output,*)' '
!  write(i_unit_output,*)'Nucl_num_shell_Aos(i) = ',Nucl_num_shell_Aos(i)
   do j = 1, Nucl_num_shell_Aos(i)
    i_shell +=1 
    i_ao = Nucl_list_shell_Aos(i,j)
    character_shell = trim(ao_l_char(i_ao))
!   write(i_unit_output,*),j,i_shell,i_ao!trim(character_shell)
    do k = 1, ao_prim_num(i_ao)
     i_prim +=1
     if(i_prim.lt.100)then
       write(i_unit_output,'(4X,I3,3X,A1,6X,I2,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo(i_ao,k),ao_coef(i_ao,k)
     else
       write(i_unit_output,'(4X,I3,3X,A1,5X,I3,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo(i_ao,k),ao_coef(i_ao,k)
     endif
    enddo
    write(i_unit_output,*)''
   enddo
 enddo

 write(i_unit_output,*)''
 write(i_unit_output,'(A47,2X,I3)')'TOTAL NUMBER OF BASIS SET SHELLS             =', i_shell
 write(i_unit_output,'(A47,2X,I3)')'NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =', ao_num
! this is for the new version of molden
 write(i_unit_output,'(A12)')'PP    =NONE'    
 write(i_unit_output,*)''


end

subroutine write_Mo_basis(i_unit_output)
 implicit none
 integer, intent(in) :: i_unit_output
 integer :: i,j,k,l, getUnitAndOpen
 integer :: i_5,i_mod

 write(i_unit_output,*),'          ----------------------'
 write(i_unit_output,*),'          MCSCF NATURAL ORBITALS'
 write(i_unit_output,*),'          ----------------------'
 write(i_unit_output,*),'                                '

 do j = 1, mo_tot_num
  write(i_unit_output,'(18X,I3)')j
  write(i_unit_output,*)''
  write(i_unit_output,'(18X,F8.5)')-1.d0
  write(i_unit_output,*)''
  do i = 1, ao_num
!   write(i_unit_output,'(2X,I3, 2X A1,  I3, 2X A4 , F9.6)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i)),mo_coef(i,j)
! F12.6 for larger coefficients...
   write(i_unit_output,'(2X,I3, 2X A1,  I3, 2X A4 , F12.6)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i)),mo_coef(i,j)
!  write(i_unit_output,'(I3, X A1, X I3, X A4 X F16.8)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i))
  enddo
  write(i_unit_output,*)''
 enddo

end