An ebuild for abinit-6.4.2 using external dependencies as proper Gentoo packages. Based on bug 249493, still work in progress. With USE="smp" compiles, but crashes.

1 files changed, 3096 insertions, 0 deletions

diff --git a/sci-physics/abinit/files/6.4.2-openmp.patch b/sci-physics/abinit/files/6.4.2-openmp.patch
new file mode 100644
index 000000000..e3e882ecf
--- /dev/null
+++ b/sci-physics/abinit/files/6.4.2-openmp.patch
@@ -0,0 +1,3096 @@
+diff -Naur abinit-6.4.2.bak/src/10_defs/defs_datatypes.F90 abinit-6.4.2/src/10_defs/defs_datatypes.F90
+--- src/10_defs/defs_datatypes.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/10_defs/defs_datatypes.F90	2011-01-05 09:39:41.000000000 +0000
+@@ -169,10 +169,10 @@
+   integer :: nsppol                ! number of spin-polarizations
+   integer :: occopt                ! Occupation option, see input variable.
+ 
+-  !$integer :: kptopt
+-  !$real(dp) :: tolwfr
+-  !$real(dp),pointer :: resid(mband*nkpt*nsppol)
+-  !$resid(mband*nkpt*nsppol)=residuals (hartree**2)
++  !%integer :: kptopt
++  !%real(dp) :: tolwfr
++  !%real(dp),pointer :: resid(mband*nkpt*nsppol)
++  !%resid(mband*nkpt*nsppol)=residuals (hartree**2)
+ 
+   real(dp) :: entropy              ! Entropy associated with the smearing (adimensional)
+   real(dp) :: fermie               ! Fermi energy
+@@ -278,7 +278,7 @@
+ 
+  type coeff2_type
+ 
+-  !$integer :: size(2)
++  !%integer :: size(2)
+   real(dp), pointer :: value(:,:)   !SET2NULL
+ 
+  end type coeff2_type
+@@ -295,7 +295,7 @@
+ 
+  type coeff2c_type
+ 
+-  !$integer :: size(2)
++  !%integer :: size(2)
+   complex(dpc), pointer :: value(:,:)  !SET2NULL
+ 
+  end type coeff2c_type
+@@ -312,7 +312,7 @@
+ 
+  type coeff3_type
+ 
+-  !$integer :: size(3)
++  !%integer :: size(3)
+   real(dp), pointer :: value(:,:,:)  !SET2NULL
+ 
+  end type coeff3_type
+@@ -329,7 +329,7 @@
+ 
+  type coeff4_type
+ 
+-  !$integer :: size(4)
++  !%integer :: size(4)
+   real(dp), pointer :: value(:,:,:,:)   !SET2NULL
+ 
+  end type coeff4_type
+@@ -348,7 +348,7 @@
+ 
+  type coeff5_type
+ 
+-  !$integer :: size(5)
++  !%integer :: size(5)
+   real(dp), pointer :: value(:,:,:,:,:)   !SET2NULL
+ 
+  end type coeff5_type
+@@ -1757,7 +1757,7 @@
+   integer :: cplex_dij
+    ! cplex=1 if dij are real, 2 if they are complex
+ 
+-  !$integer :: has_dijexxcore
++  !%integer :: has_dijexxcore
+    ! 1 if dijexxcore is allocated
+    ! 2 if dijexxcore is already computed
+ 
+@@ -1839,7 +1839,7 @@
+    !  dij(:,:,3) contains Dij^up-dn (only if nspinor=2)
+    !  dij(:,:,4) contains Dij^dn-up (only if nspinor=2)
+ 
+-  !$real(dp),pointer :: dijexxcore(:,:)
++  !%real(dp),pointer :: dijexxcore(:,:)
+   ! dijexxcore(cplex_dij*lmn2_size,ndij)
+   ! Onsite matrix elements of the Fock operator generated by core electrons
+ 
+@@ -2551,21 +2551,21 @@
+ 
+ ! All the energies are in Hartree, obtained "per unit cell".
+   type(energies_type) :: energies
+-!!$  real(dp) :: eei      ! local pseudopotential energy (Hartree)
+-!!$  real(dp) :: eeig     ! sum of eigenvalue energy (Hartree)
+-!!$  real(dp) :: eew      ! Ewald energy (Hartree)
+-!!$  real(dp) :: ehart    ! Hartree part of total energy (Hartree)
+-!!$  real(dp) :: eii      ! pseudopotential core-core energy
+-!!$  real(dp) :: ek       ! kinetic energy (Hartree)
+-!!$  real(dp) :: enefield ! the term of the energy functional that depends
+-!!$                       ! explicitely on the electric field
+-!!$                       ! enefield = -ucvol*E*P
+-!!$  real(dp) :: enl      ! nonlocal pseudopotential energy (Hartree)
++!!%  real(dp) :: eei      ! local pseudopotential energy (Hartree)
++!!%  real(dp) :: eeig     ! sum of eigenvalue energy (Hartree)
++!!%  real(dp) :: eew      ! Ewald energy (Hartree)
++!!%  real(dp) :: ehart    ! Hartree part of total energy (Hartree)
++!!%  real(dp) :: eii      ! pseudopotential core-core energy
++!!%  real(dp) :: ek       ! kinetic energy (Hartree)
++!!%  real(dp) :: enefield ! the term of the energy functional that depends
++!!%                       ! explicitely on the electric field
++!!%                       ! enefield = -ucvol*E*P
++!!%  real(dp) :: enl      ! nonlocal pseudopotential energy (Hartree)
+   real(dp) :: entropy  ! entropy (Hartree)
+-!!$  real(dp) :: enxc     ! exchange-correlation energy (Hartree)
+-!!$  real(dp) :: enxcdc   ! exchange-correlation double-counting energy (Hartree)
+-!!$  real(dp) :: epaw     ! PAW spherical energy (Hartree)
+-!!$  real(dp) :: epawdc   ! PAW spherical double-counting energy (Hartree)
++!!%  real(dp) :: enxc     ! exchange-correlation energy (Hartree)
++!!%  real(dp) :: enxcdc   ! exchange-correlation double-counting energy (Hartree)
++!!%  real(dp) :: epaw     ! PAW spherical energy (Hartree)
++!!%  real(dp) :: epawdc   ! PAW spherical double-counting energy (Hartree)
+   real(dp) :: etotal   ! total energy (Hartree)
+                        ! for fixed occupation numbers (occopt==0,1,or 2):
+                        !   etotal=ek+ehart+enxc+eei+eew+eii+enl+PAW_spherical_part
+@@ -3009,7 +3009,7 @@
+ ! WARNING : if you modify this datatype, please check there there is no creation/destruction/copy routine,
+ ! declared in another part of ABINIT, that might need to take into account your modification.
+ 
+-  integer :: stat = 1 !$=JPT_ISPOINTER  =1
++  integer :: stat = 1 !%=JPT_ISPOINTER  =1
+   complex(gwpc),pointer :: datum(:,:)  SET2NULL
+  end type Jpt_gwpc_2D
+ !!***
+@@ -3034,7 +3034,7 @@
+ ! WARNING : if you modify this datatype, please check there there is no creation/destruction routine,
+ ! declared in another part of ABINIT, that might need to take into account your modification.
+ 
+-  integer :: stat = 1 !$=JPT_ISPOINTER  =1
++  integer :: stat = 1 !%=JPT_ISPOINTER  =1
+   complex(gwpc),pointer :: datum(:,:,:)  SET2NULL
+  end type Jpt_gwpc_3D
+ !!***
+@@ -3071,8 +3071,8 @@
+    logical :: save_memory_devel=.FALSE.
+ 
+    ! Use in-place storage of the PPm parameters, should be FALSE for non-developers.
+-   !$real(dp) :: zcut
+-   !$real(dp),pointer :: qibz(:,:)
++   !%real(dp) :: zcut
++   !%real(dp),pointer :: qibz(:,:)
+ 
+    !logical :: has_inversion
+    !logical :: has_time_reversal
+diff -Naur abinit-6.4.2.bak/src/10_defs/defs_fftdata.F90 abinit-6.4.2/src/10_defs/defs_fftdata.F90
+--- src/10_defs/defs_fftdata.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/10_defs/defs_fftdata.F90	2011-01-05 09:40:20.000000000 +0000
+@@ -181,8 +181,8 @@
+   write(msg,'(4a,i8)')ch10,&
+ &  ' size_goed_fft : BUG-',ch10,&
+ &  ' nbest = ',nbest 
+-  !$call wrtout(std_out,msg,'COLL') 
+-  !$call leave_new('COLL')
++  !%call wrtout(std_out,msg,'COLL') 
++  !%call leave_new('COLL')
+   write(std_out,*)msg
+   STOP
+  end if
+@@ -192,8 +192,8 @@
+ &  ' size_goed_fft : ERROR-',ch10,&
+ &  ' nbest = ',nbest,ch10,&
+ &  ' is larger than any allowable FFT'
+-  !$call wrtout(std_out,msg,'COLL') 
+-  !$call leave_new('COLL')
++  !%call wrtout(std_out,msg,'COLL') 
++  !%call leave_new('COLL')
+   write(std_out,*)msg
+   STOP
+  end if
+diff -Naur abinit-6.4.2.bak/src/12_hide_mpi/m_xmpi.F90 abinit-6.4.2/src/12_hide_mpi/m_xmpi.F90
+--- src/12_hide_mpi/m_xmpi.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/12_hide_mpi/m_xmpi.F90	2011-01-05 09:02:11.000000000 +0000
+@@ -86,7 +86,7 @@
+   integer :: nprocs = xmpi_undefined
+   ! The number of processors in the communicator.
+ 
+-  !$integer,pointer :: ranks_in_world(:)   SET2NULL
++  !%integer,pointer :: ranks_in_world(:)   SET2NULL
+   ! The MPI ranks in MPI_COMM_WORLD of the nodes beloging to the communicator.
+ 
+  end type mpicomm_t
+diff -Naur abinit-6.4.2.bak/src/18_timing/m_timer.F90 abinit-6.4.2/src/18_timing/m_timer.F90
+--- src/18_timing/m_timer.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/18_timing/m_timer.F90	2011-01-05 17:23:06.000000000 +0000
+@@ -21,7 +21,7 @@
+ MODULE m_timer
+ 
+  use defs_basis
+- !$use m_errors TODO move m_timer to a lower level
++ !%use m_errors TODO move m_timer to a lower level
+ 
+  implicit none
+ 
+@@ -261,7 +261,7 @@
+  CASE DEFAULT
+    write(std_out,'(a,i0)')" Wrong value for topt: ",topt
+    stop
+-   !$MSG_ERROR(msg)
++   !%MSG_ERROR(msg)
+  END SELECT
+ 
+ end subroutine timing
+diff -Naur abinit-6.4.2.bak/src/27_toolbox_oop/m_errors.F90 abinit-6.4.2/src/27_toolbox_oop/m_errors.F90
+--- src/27_toolbox_oop/m_errors.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/27_toolbox_oop/m_errors.F90	2011-01-05 17:33:42.000000000 +0000
+@@ -748,7 +748,7 @@
+ &  TRIM(my_msg),toupper(level),ch10,&
+ &  TRIM(message)
+   call wrtout(std_out,sbuf,mode_paral) 
+-  !$call dump_config(std_out)
++  !%call dump_config(std_out)
+   call dump_config_fake()
+   call leave_new('COLL')
+ 
+diff -Naur abinit-6.4.2.bak/src/32_util/m_radmesh.F90 abinit-6.4.2/src/32_util/m_radmesh.F90
+--- src/32_util/m_radmesh.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/32_util/m_radmesh.F90	2011-01-05 08:51:17.000000000 +0000
+@@ -451,7 +451,7 @@
+  !gg = zero
+ 
+  int_meshsz=Pawrad%int_meshsz
+- !$int_meshsz=Pawrad%mesh_size
++ !%int_meshsz=Pawrad%mesh_size
+ 
+  ! the line below requires hh as work array. 
+  hh = ff2  
+diff -Naur abinit-6.4.2.bak/src/32_util/m_special_funcs.F90 abinit-6.4.2/src/32_util/m_special_funcs.F90
+--- src/32_util/m_special_funcs.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/32_util/m_special_funcs.F90	2011-01-05 08:51:59.000000000 +0000
+@@ -232,7 +232,7 @@
+ 
+  r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
+  if (r<PPAD) r=r+PPAD
+- !$if (r<tol10) RETURN
++ !%if (r<tol10) RETURN
+ 
+  rxy=SQRT(kcart(1)**2+kcart(2)**2)
+  if (rxy<PPAD)rxy=r+PPAD
+@@ -391,7 +391,7 @@
+ 
+  r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
+  if (r<PPAD) r=r+PPAD
+- !$if (r<tol10) RETURN
++ !%if (r<tol10) RETURN
+ 
+  rxy=SQRT(kcart(1)**2+kcart(2)**2)
+  if (rxy<PPAD) rxy=r+PPAD
+diff -Naur abinit-6.4.2.bak/src/42_geometry/m_crystal.F90 abinit-6.4.2/src/42_geometry/m_crystal.F90
+--- src/42_geometry/m_crystal.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/42_geometry/m_crystal.F90	2011-01-05 09:25:18.000000000 +0000
+@@ -76,7 +76,7 @@
+   integer :: ntypat   
+   ! Number of type of atoms
+ 
+-  !$integer :: ntypalch,ntyppure
++  !%integer :: ntypalch,ntyppure
+ 
+   integer :: npsp  
+   ! No. of pseudopotentials
+diff -Naur abinit-6.4.2.bak/src/43_ptgroups/m_defs_ptgroups.F90 abinit-6.4.2/src/43_ptgroups/m_defs_ptgroups.F90
+--- src/43_ptgroups/m_defs_ptgroups.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/43_ptgroups/m_defs_ptgroups.F90	2011-01-05 16:44:37.000000000 +0000
+@@ -147,11 +147,11 @@
+    ! NB: operations are referred to the standard coordinate system.
+    ! Page 815-816 of International Tables for crystallography Vol.A.
+ 
+-   !$integer,pointer :: symafm(:)   SET2NULL
++   !%integer,pointer :: symafm(:)   SET2NULL
+    ! symafm(nsym)
+    ! AFM part of the symmetry operation
+ 
+-   !$real(dp),pointer :: tnons(:,:)  SET2NULL
++   !%real(dp),pointer :: tnons(:,:)  SET2NULL
+    ! tnons(3,nsym)
+    ! fractional translations.
+ 
+diff -Naur abinit-6.4.2.bak/src/43_ptgroups/m_ptgroups.F90 abinit-6.4.2/src/43_ptgroups/m_ptgroups.F90
+--- src/43_ptgroups/m_ptgroups.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/43_ptgroups/m_ptgroups.F90	2011-01-05 16:45:10.000000000 +0000
+@@ -894,7 +894,7 @@
+    read(unt,*,ERR=10) nirreps_k
+    ABI_CHECK(Gk%nclass == nirreps_k,"Gk%nclass /= nirreps_k")
+ 
+-   !$$ allocate(Gk%class_names(Gk%nclass)) 
++   !%$$ allocate(Gk%class_names(Gk%nclass)) 
+    allocate(Gk%Irreps(nirreps_k))
+ 
+    do irp=1,nirreps_k
+diff -Naur abinit-6.4.2.bak/src/49_gw_toolbox_oop/m_bs_defs.F90 abinit-6.4.2/src/49_gw_toolbox_oop/m_bs_defs.F90
+--- src/49_gw_toolbox_oop/m_bs_defs.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/49_gw_toolbox_oop/m_bs_defs.F90	2011-01-05 16:37:04.000000000 +0000
+@@ -43,23 +43,23 @@
+ !!***
+ 
+ ! Approximation for the excitonic Hamiltonian
+-!$ integer,public,parameter :: BSE_HRPA = 1   
+-!$ integer,public,parameter :: BSE_HGW  = 2
+-!$ integer,public,parameter :: BSE_HEXC = 3
++!% integer,public,parameter :: BSE_HRPA = 1   
++!% integer,public,parameter :: BSE_HGW  = 2
++!% integer,public,parameter :: BSE_HEXC = 3
+ 
+ ! Treatment of the resontant and coupling blocks of the Hamiltonian
+-!$? integer,public,parameter :: BSE_HRPA = 1   
+-!$? integer,public,parameter :: BSE_HGW  = 2
++!%? integer,public,parameter :: BSE_HRPA = 1   
++!%? integer,public,parameter :: BSE_HGW  = 2
+ 
+-!$
+-!$! Treatment of W(G1,G2)
+-!$ integer,public,parameter :: BSE_WGG_DIAGONAL=1
+-!$ integer,public,parameter :: BSE_WGG_FULL    =2
+-
+-!$! Treatment of W(omega)
+-!$ integer,public,parameter :: BSE_WFREQ_STATIC=1
+-!$ integer,public,parameter :: BSE_WFREQ_PPM   =2
+-!$ integer,public,parameter :: BSE_WFREQ_FULL  =2
++!%
++!%! Treatment of W(G1,G2)
++!% integer,public,parameter :: BSE_WGG_DIAGONAL=1
++!% integer,public,parameter :: BSE_WGG_FULL    =2
++
++!%! Treatment of W(omega)
++!% integer,public,parameter :: BSE_WFREQ_STATIC=1
++!% integer,public,parameter :: BSE_WFREQ_PPM   =2
++!% integer,public,parameter :: BSE_WFREQ_FULL  =2
+ 
+ !----------------------------------------------------------------------
+ 
+@@ -85,7 +85,7 @@
+   logical :: wdiag             ! Use diagonal approximation for W.
+   logical :: wfull             ! Use full inverse dielectric matrix.
+   logical :: use_haydock       ! Use haydock iterative method to calculate epsilon.
+-  !$logical :: use_cg          ! Use conjugate gradient minimization techniques to calculate the eigenstates.
++  !%$logical :: use_cg          ! Use conjugate gradient minimization techniques to calculate the eigenstates.
+   logical :: gwterm            ! Add QP corrections to the band structure either with scissors or GW file.
+   logical :: exchangeterm      ! Include exchange terms in the BS Hamiltonian.
+   logical :: coulombterm       ! Include W term in the BS Hamiltonian.
+@@ -99,7 +99,7 @@
+ 
+   integer :: npweps            ! No. of G in the Screening.
+   integer :: npwwfn            ! No. of G for wave functions.
+-  !$integer :: npwx              ! No. of G for the exchange part.
++  !%integer :: npwx              ! No. of G for the exchange part.
+   integer :: npwvec            ! MAX between npwwfn and npweps
+   integer :: nsh               ! Number of shells corresponding to npwvec=MAX(npwwfn,npweps)
+   integer :: nbnds             ! Total number of bands considered.
+diff -Naur abinit-6.4.2.bak/src/51_manage_mpi/mpi_enreg_tools.F90 abinit-6.4.2/src/51_manage_mpi/mpi_enreg_tools.F90
+--- src/51_manage_mpi/mpi_enreg_tools.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/51_manage_mpi/mpi_enreg_tools.F90	2011-01-05 08:53:13.000000000 +0000
+@@ -675,8 +675,8 @@
+    MSG_ERROR(" %paral_compil_kpt/=1")
+  end if
+ 
+-!$spaceComm=MPI_enreg%spaceComm
+-!$nprocs = xcomm_size(spaceComm)
++!%spaceComm=MPI_enreg%spaceComm
++!%nprocs = xcomm_size(spaceComm)
+  nprocs  = MPI_enreg%nproc
+  my_rank = MPI_enreg%me
+ 
+diff -Naur abinit-6.4.2.bak/src/52_fft_mpi_noabirule/accrho.F90 abinit-6.4.2/src/52_fft_mpi_noabirule/accrho.F90
+--- src/52_fft_mpi_noabirule/accrho.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/52_fft_mpi_noabirule/accrho.F90	2011-01-06 16:08:27.000000000 +0000
+@@ -91,14 +91,14 @@
+         integer unused
+         unused=0
+ ! *************************************************************************
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_num_threads ( )
+-!$        end function omp_get_num_threads
+-!$      end interface
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_thread_num ( )
+-!$        end function omp_get_thread_num
+-!$      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_num_threads ( )
++!%        end function omp_get_num_threads
++!%      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_thread_num ( )
++!%        end function omp_get_thread_num
++!%      end interface
+ 
+         write(6,*)' accrho : enter '
+ 
+diff -Naur abinit-6.4.2.bak/src/52_fft_mpi_noabirule/applypot.F90 abinit-6.4.2/src/52_fft_mpi_noabirule/applypot.F90
+--- src/52_fft_mpi_noabirule/applypot.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/52_fft_mpi_noabirule/applypot.F90	2011-01-06 16:09:25.000000000 +0000
+@@ -91,14 +91,14 @@
+     integer unused
+         unused=0
+ 
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_num_threads ( )
+-!$        end function omp_get_num_threads
+-!$      end interface
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_thread_num ( )
+-!$        end function omp_get_thread_num
+-!$      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_num_threads ( )
++!%        end function omp_get_num_threads
++!%      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_thread_num ( )
++!%        end function omp_get_thread_num
++!%      end interface
+ 
+         write(6,*)' applypot : enter '
+ 
+diff -Naur abinit-6.4.2.bak/src/52_fft_mpi_noabirule/back.F90 abinit-6.4.2/src/52_fft_mpi_noabirule/back.F90
+--- src/52_fft_mpi_noabirule/back.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/52_fft_mpi_noabirule/back.F90	2011-01-07 14:21:34.000000000 +0000
+@@ -90,14 +90,14 @@
+ #endif
+ ! *************************************************************************
+ 
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_num_threads ( )
+-!$        end function omp_get_num_threads
+-!$      end interface
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_thread_num ( )
+-!$        end function omp_get_thread_num
+-!$      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_num_threads ( )
++!%        end function omp_get_num_threads
++!%      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_thread_num ( )
++!%        end function omp_get_thread_num
++!%      end interface
+ 
+ !DEBUG
+ !       write(6,*)' back : enter '
+@@ -121,7 +121,8 @@
+ 
+         lock=0
+ !$omp parallel  default(private) &
+-!$omp shared(ndat,n1,n2,n3,nd1,nd2,nd3,nd2proc,nd3proc,mpi_enreg%me_fft,mpi_enreg%nproc_fft,ncache,zr,zf,lock,icplex)
++!%omp shared(ndat,n1,n2,n3,nd1,nd2,nd3,nd2proc,nd3proc,mpi_enreg%me_fft,mpi_enreg%nproc_fft,ncache,zr,zf,lock,icplex)
++!$omp shared(ndat,n1,n2,n3,nd1,nd2,nd3,nd2proc,nd3proc,mpi_enreg,ncache,zr,zf,lock,icplex)
+ 
+         iam=0
+         npr=1
+diff -Naur abinit-6.4.2.bak/src/52_fft_mpi_noabirule/back_wf.F90 abinit-6.4.2/src/52_fft_mpi_noabirule/back_wf.F90
+--- src/52_fft_mpi_noabirule/back_wf.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/52_fft_mpi_noabirule/back_wf.F90	2011-01-07 14:18:09.000000000 +0000
+@@ -98,14 +98,14 @@
+         real(KIND=8) :: tsec(2)
+         type(MPI_type),intent(inout) :: mpi_enreg
+         integer :: old_paral_level
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_num_threads ( )
+-!$        end function omp_get_num_threads
+-!$      end interface
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_thread_num ( )
+-!$        end function omp_get_thread_num
+-!$      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_num_threads ( )
++!%        end function omp_get_num_threads
++!%      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_thread_num ( )
++!%        end function omp_get_thread_num
++!%      end interface
+ 
+         if(.false.)write(6,*)iproc,nproc
+         if(mpi_enreg%mode_para=='b') ioption=1
+diff -Naur abinit-6.4.2.bak/src/52_fft_mpi_noabirule/forw.F90 abinit-6.4.2/src/52_fft_mpi_noabirule/forw.F90
+--- src/52_fft_mpi_noabirule/forw.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/52_fft_mpi_noabirule/forw.F90	2011-01-07 14:21:49.000000000 +0000
+@@ -95,14 +95,14 @@
+         integer :: old_paral_level
+ #endif
+ ! *************************************************************************
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_num_threads ( )
+-!$        end function omp_get_num_threads
+-!$      end interface
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_thread_num ( )
+-!$        end function omp_get_thread_num
+-!$      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_num_threads ( )
++!%        end function omp_get_num_threads
++!%      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_thread_num ( )
++!%        end function omp_get_thread_num
++!%      end interface
+ 
+ !DEBUG
+ !       write(6,*)' forw : enter '
+@@ -134,7 +134,8 @@
+ 
+         lock=0
+ !$omp parallel  default(private) &
+-!$omp shared(ndat,n1,n2,n3,nd1,nd2,nd3,nd2proc,nd3proc,mpi_enreg%me_fft,mpi_enreg%nproc_fft,ncache,zr,zf,lock,icplex)
++!%omp shared(ndat,n1,n2,n3,nd1,nd2,nd3,nd2proc,nd3proc,mpi_enreg%me_fft,mpi_enreg%nproc_fft,ncache,zr,zf,lock,icplex)
++!$omp shared(ndat,n1,n2,n3,nd1,nd2,nd3,nd2proc,nd3proc,mpi_enreg,ncache,zr,zf,lock,icplex)
+ 
+         iam=0
+         npr=1
+diff -Naur abinit-6.4.2.bak/src/52_fft_mpi_noabirule/forw_wf.F90 abinit-6.4.2/src/52_fft_mpi_noabirule/forw_wf.F90
+--- src/52_fft_mpi_noabirule/forw_wf.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/52_fft_mpi_noabirule/forw_wf.F90	2011-01-07 14:18:58.000000000 +0000
+@@ -98,14 +98,14 @@
+         type(MPI_type),intent(inout) :: mpi_enreg
+          integer :: old_paral_level
+ 
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_num_threads ( )
+-!$        end function omp_get_num_threads
+-!$      end interface
+-!$      interface
+-!$        integer ( kind=4 ) function omp_get_thread_num ( )
+-!$        end function omp_get_thread_num
+-!$      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_num_threads ( )
++!%        end function omp_get_num_threads
++!%      end interface
++!%      interface
++!%        integer ( kind=4 ) function omp_get_thread_num ( )
++!%        end function omp_get_thread_num
++!%      end interface
+ 
+         if(mpi_enreg%mode_para=='b') ioption=1
+ !       call timab(542,1,tsec)
+diff -Naur abinit-6.4.2.bak/src/52_fft_mpi_noabirule/m_fftw3.F90 abinit-6.4.2/src/52_fft_mpi_noabirule/m_fftw3.F90
+--- src/52_fft_mpi_noabirule/m_fftw3.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/52_fft_mpi_noabirule/m_fftw3.F90	2011-01-05 09:58:51.000000000 +0000
+@@ -156,7 +156,7 @@
+    !integer,pointer :: n(:) rank
+    !integer,pointer :: inembed(:) rank     
+    !integer,pointer :: onembed(:) rank
+-   !$integer(C_INT) :: alignment(2)                ! The alignment of the arrays used to costruct the plan.
++   !%integer(C_INT) :: alignment(2)                ! The alignment of the arrays used to costruct the plan.
+                                                    ! needed because we cannot use fftw_malloc.
+    integer :: ndiv(3)=-1                           ! The number of FFT divisions.
+  end type fftw_plan_type
+@@ -1022,10 +1022,10 @@
+  allocate(ff(nx*ny*nz*ndat),gg(2,nx*ny*nz*ndat))
+ 
+  ! These calls seem to make the code stuck if we link against MKL
+- !$call fftw3_r2c_op(nx,ny,nz,nx,ny,nz,ndat,ff,gg,fftw_flags=wisdom_flag)
+- !$call fftw3_c2r_op(nx,ny,nz,nx,ny,nz,ndat,ff,gg,fftw_flags=wisdom_flag)
+- !$write(msg,'(a)')" fftw3_c2r_op done"
+- !$call wrtout(std_out,msg,"COLL")
++ !%call fftw3_r2c_op(nx,ny,nz,nx,ny,nz,ndat,ff,gg,fftw_flags=wisdom_flag)
++ !%call fftw3_c2r_op(nx,ny,nz,nx,ny,nz,ndat,ff,gg,fftw_flags=wisdom_flag)
++ !%write(msg,'(a)')" fftw3_c2r_op done"
++ !%call wrtout(std_out,msg,"COLL")
+ 
+  deallocate(ff,gg)
+ 
+diff -Naur abinit-6.4.2.bak/src/53_abiutil/m_abilasi.F90 abinit-6.4.2/src/53_abiutil/m_abilasi.F90
+--- src/53_abiutil/m_abilasi.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/53_abiutil/m_abilasi.F90	2011-01-05 17:36:48.000000000 +0000
+@@ -567,36 +567,36 @@
+ 
+   MSG_PERS_ERROR("Not coded yet")
+ 
+-!$$  call init_scalapack(Slk_processor,comm)
+-!$$  istwf_k=1
+-!$$
+-!$$  ! Initialize and fill Scalapack matrix from the global one.
+-!$$  tbloc=SLK_BLOCK_SIZE
+-!$$  call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
+-!$$
+-!$$  write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
+-!$$  call wrtout(std_out,msg,"PERS")
+-!$$
+-!$$  call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
+-!$$
+-!$$  want_eigenvectors = starts_with(jobz,(/"V","v"/))
+-!$$  if (want_eigenvectors) then ! Initialize the distributed vectors.
+-!$$   call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
+-!$$  end if
+-!$$
+-!$$  ! Solve the problem with scaLAPACK.
+-!$$  call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w)
+-!$$
+-!$$  call destruction_matrix_scalapack(Slk_mat)
+-!$$  
+-!$$  if (want_eigenvectors) then ! A is overwritten with the eigenvectors 
+-!$$   a = czero  
+-!$$   call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node.
+-!$$   call destruction_matrix_scalapack(Slk_vec)
+-!$$   call xsum_mpi(a,comm,ierr)                        ! Fill the remaing entries of the global matrix
+-!$$  end if
+-!$$
+-!$$  call end_scalapack(Slk_processor)
++!%%  call init_scalapack(Slk_processor,comm)
++!%%  istwf_k=1
++!%%
++!%%  ! Initialize and fill Scalapack matrix from the global one.
++!%%  tbloc=SLK_BLOCK_SIZE
++!%%  call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
++!%%
++!%%  write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
++!%%  call wrtout(std_out,msg,"PERS")
++!%%
++!%%  call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
++!%%
++!%%  want_eigenvectors = starts_with(jobz,(/"V","v"/))
++!%%  if (want_eigenvectors) then ! Initialize the distributed vectors.
++!%%   call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
++!%%  end if
++!%%
++!%%  ! Solve the problem with scaLAPACK.
++!%%  call slk_pzheev(jobz,uplo,Slk_mat,Slk_vec,w)
++!%%
++!%%  call destruction_matrix_scalapack(Slk_mat)
++!%%  
++!%%  if (want_eigenvectors) then ! A is overwritten with the eigenvectors 
++!%%   a = czero  
++!%%   call slk_matrix_to_global_dpc_2D(Slk_vec,"All",a) ! Fill the entries calculated by this node.
++!%%   call destruction_matrix_scalapack(Slk_vec)
++!%%   call xsum_mpi(a,comm,ierr)                        ! Fill the remaing entries of the global matrix
++!%%  end if
++!%%
++!%%  call end_scalapack(Slk_processor)
+ 
+   RETURN 
+ #endif
+@@ -1086,7 +1086,7 @@
+   ! Solve the problem with scaLAPACK.
+   MSG_ERROR("slk_pZHEGV not yet coded")
+   ! TODO
+-  !$call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
++  !%call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
+ 
+   call destruction_matrix_scalapack(Slk_matB)
+   
+@@ -1324,37 +1324,37 @@
+ 
+   MSG_ERROR("Not coded yet")
+ 
+-  !$$ call init_scalapack(Slk_processor,comm)
+-  !$$ istwf_k=1
++  !%% call init_scalapack(Slk_processor,comm)
++  !%% istwf_k=1
+ 
+-  !$$ ! Initialize and fill Scalapack matrix from the global one.
+-  !$$ tbloc=SLK_BLOCK_SIZE
++  !%% ! Initialize and fill Scalapack matrix from the global one.
++  !%% tbloc=SLK_BLOCK_SIZE
+ 
+-  !$$ write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
+-  !$$ call wrtout(std_out,msg,"PERS")
++  !%% write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
++  !%% call wrtout(std_out,msg,"PERS")
+ 
+-  !$$ call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
+-  !$$ call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
++  !%% call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
++  !%% call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
+ 
+-  !$$ call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
+-  !$$ call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
++  !%% call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
++  !%% call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
+ 
+-  !$$ ! Solve the problem with scaLAPACK.
+-  !$$ MSG_ERROR("slk_pZHEGV not yet coded")
+-  !$$ ! TODO
+-  !$$ !$call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
++  !%% ! Solve the problem with scaLAPACK.
++  !%% MSG_ERROR("slk_pZHEGV not yet coded")
++  !%% ! TODO
++  !%% !%call slk_pzhegv(itype,jobz,uplo,Slk_matA,Slk_matB,w)
+ 
+-  !$$ call destruction_matrix_scalapack(Slk_matB)
+-  !$$ 
+-  !$$ if (starts_with(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors 
+-  !$$  a = czero  
+-  !$$  call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node.
+-  !$$  call xsum_mpi(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
+-  !$$ end if
++  !%% call destruction_matrix_scalapack(Slk_matB)
++  !%% 
++  !%% if (starts_with(jobz,(/"V","v"/))) then ! A is overwritten with the eigenvectors 
++  !%%  a = czero  
++  !%%  call slk_matrix_to_global_dpc_2D(Slk_matA,"All",a) ! Fill the entries calculated by this node.
++  !%%  call xsum_mpi(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
++  !%% end if
+ 
+-  !$$ call destruction_matrix_scalapack(Slk_matA)
++  !%% call destruction_matrix_scalapack(Slk_matA)
+ 
+-  !$$ call end_scalapack(Slk_processor)
++  !%% call end_scalapack(Slk_processor)
+ 
+   RETURN 
+ #endif
+@@ -1868,36 +1868,36 @@
+ #if defined HAVE_LINALG_MPI
+ 
+   MSG_ERROR("Not coded yet")
+-  !$$ call init_scalapack(Slk_processor,comm)
+-  !$$ istwf_k=1
++  !%% call init_scalapack(Slk_processor,comm)
++  !%% istwf_k=1
+ 
+-  !$$ ! Initialize and fill Scalapack matrix from the global one.
+-  !$$ tbloc=SLK_BLOCK_SIZE
+-  !$$ call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
+-
+-  !$$ write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
+-  !$$ call wrtout(std_out,msg,"PERS")
+-
+-  !$$ call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
+-
+-  !$$ want_eigenvectors = starts_with(jobz,(/"V","v"/))
+-  !$$ if (want_eigenvectors) then ! Initialize the distributed vectors.
+-  !$$  call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
+-  !$$ end if
+-
+-  !$$ ! Solve the problem.
+-  !$$ call slk_pzheevx(jobz,range,uplo,Slk_mat,vl,vu,il,iu,abstol,Slk_vec,m,w)
+-
+-  !$$ call destruction_matrix_scalapack(Slk_mat)
+-  !$$ 
+-  !$$ if (want_eigenvectors) then ! A is overwritten with the eigenvectors 
+-  !$$  z = czero  
+-  !$$  call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
+-  !$$  call destruction_matrix_scalapack(Slk_vec)
+-  !$$  call xsum_mpi(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
+-  !$$ end if
++  !%% ! Initialize and fill Scalapack matrix from the global one.
++  !%% tbloc=SLK_BLOCK_SIZE
++  !%% call init_matrix_scalapack(Slk_mat,n,n,Slk_processor,istwf_k,tbloc=tbloc)
++
++  !%% write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
++  !%% call wrtout(std_out,msg,"PERS")
++
++  !%% call slk_matrix_from_global_dpc_2D(Slk_mat,uplo,a)
++
++  !%% want_eigenvectors = starts_with(jobz,(/"V","v"/))
++  !%% if (want_eigenvectors) then ! Initialize the distributed vectors.
++  !%%  call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
++  !%% end if
++
++  !%% ! Solve the problem.
++  !%% call slk_pzheevx(jobz,range,uplo,Slk_mat,vl,vu,il,iu,abstol,Slk_vec,m,w)
++
++  !%% call destruction_matrix_scalapack(Slk_mat)
++  !%% 
++  !%% if (want_eigenvectors) then ! A is overwritten with the eigenvectors 
++  !%%  z = czero  
++  !%%  call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
++  !%%  call destruction_matrix_scalapack(Slk_vec)
++  !%%  call xsum_mpi(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
++  !%% end if
+ 
+-  !$$ call end_scalapack(Slk_processor)
++  !%% call end_scalapack(Slk_processor)
+ 
+   RETURN 
+ #endif
+@@ -2159,7 +2159,7 @@
+   ! Solve the problem.
+   MSG_ERROR("slk_pZHEGVX not coded yet")
+   ! TODO write the scaLAPACK wrapper.
+-  !$call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
++  !%call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
+ 
+   call destruction_matrix_scalapack(Slk_matA)
+   call destruction_matrix_scalapack(Slk_matB)
+@@ -2466,44 +2466,44 @@
+ #if defined HAVE_LINALG_MPI
+ 
+   MSG_ERROR("not coded yet")
+-  !$$ call init_scalapack(Slk_processor,comm)
+-  !$$ istwf_k=1
++  !%% call init_scalapack(Slk_processor,comm)
++  !%% istwf_k=1
+ 
+-  !$$ ! Initialize and fill Scalapack matrix from the global one.
+-  !$$ tbloc=SLK_BLOCK_SIZE
++  !%% ! Initialize and fill Scalapack matrix from the global one.
++  !%% tbloc=SLK_BLOCK_SIZE
+ 
+-  !$$ write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
+-  !$$ call wrtout(std_out,msg,"PERS")
++  !%% write(msg,'(2(a,i0))')" Using scaLAPACK version with nprocs= ",nprocs,"; block size= ",tbloc
++  !%% call wrtout(std_out,msg,"PERS")
+ 
+-  !$$ call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
+-  !$$ call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
+-
+-  !$$ call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc)
+-  !$$ call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
+-
+-  !$$ want_eigenvectors = starts_with(jobz,(/"V","v"/))
+-  !$$ if (want_eigenvectors) then ! Initialize the distributed vectors.
+-  !$$  call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
+-  !$$ end if
+-
+-  !$$ ! Solve the problem.
+-  !$$ MSG_ERROR("slk_pZHEGVX not coded yet")
+-  !$$ ! TODO write the scaLAPACK wrapper.
+-  !$$ !$call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
+-
+-  !$$ call destruction_matrix_scalapack(Slk_matA)
+-  !$$ call destruction_matrix_scalapack(Slk_matB)
+-  !$$ 
+-  !$$ if (want_eigenvectors) then ! A is overwritten with the eigenvectors 
+-  !$$  z = czero  
+-  !$$  call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
+-  !$$  call destruction_matrix_scalapack(Slk_vec)
+-  !$$  call xsum_mpi(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
+-  !$$ end if
++  !%% call init_matrix_scalapack(Slk_matA,n,n,Slk_processor,istwf_k,tbloc=tbloc)
++  !%% call slk_matrix_from_global_dpc_2D(Slk_matA,uplo,a)
++
++  !%% call init_matrix_scalapack(Slk_matB,n,n,Slk_processor,istwf_k,tbloc=tbloc)
++  !%% call slk_matrix_from_global_dpc_2D(Slk_matB,uplo,b)
++
++  !%% want_eigenvectors = starts_with(jobz,(/"V","v"/))
++  !%% if (want_eigenvectors) then ! Initialize the distributed vectors.
++  !%%  call init_matrix_scalapack(Slk_vec,n,n,Slk_processor,istwf_k,tbloc=tbloc) 
++  !%% end if
++
++  !%% ! Solve the problem.
++  !%% MSG_ERROR("slk_pZHEGVX not coded yet")
++  !%% ! TODO write the scaLAPACK wrapper.
++  !%% !%call slk_pZHEGVX(itype,jobz,range,uplo,Slk_matA,Slk_matB,vl,vu,il,iu,abstol,Slk_vec,m,w)
++
++  !%% call destruction_matrix_scalapack(Slk_matA)
++  !%% call destruction_matrix_scalapack(Slk_matB)
++  !%% 
++  !%% if (want_eigenvectors) then ! A is overwritten with the eigenvectors 
++  !%%  z = czero  
++  !%%  call slk_matrix_to_global_dpc_2D(Slk_vec,"All",z) ! Fill the entries calculated by this node.
++  !%%  call destruction_matrix_scalapack(Slk_vec)
++  !%%  call xsum_mpi(z,comm,ierr)                        ! Fill the remaing entries of the global matrix
++  !%% end if
+ 
+-  !$$ call end_scalapack(Slk_processor)
++  !%% call end_scalapack(Slk_processor)
+ 
+-  !$$ RETURN 
++  !%% RETURN 
+ #endif
+ 
+   MSG_PERS_BUG("You should not be here!")
+@@ -2927,7 +2927,7 @@
+    MSG_PERS_ERROR(msg)
+   end if
+ 
+-  !$call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a)
++  !%call slk_matrix_from_global_dpc_2D(Slk_mat,"All",a)
+ 
+   ipiv_size = my_locr(Slk_mat) + Slk_mat%descript%tab(MB_)  
+   allocate(ipiv(ipiv_size))
+@@ -2968,7 +2968,7 @@
+ 
+   ! Reconstruct the global matrix from the distributed one.
+   a = czero
+-  !$call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a)  ! Fill the entries calculated by this node.
++  !%call slk_matrix_to_global_dpc_2D(Slk_mat,"All",a)  ! Fill the entries calculated by this node.
+   call destruction_matrix_scalapack(Slk_mat)
+ 
+   call xsum_mpi(a,comm,ierr)                         ! Fill the remaing entries of the global matrix
+diff -Naur abinit-6.4.2.bak/src/53_ffts/fftw3_fourwf.F90 abinit-6.4.2/src/53_ffts/fftw3_fourwf.F90
+--- src/53_ffts/fftw3_fourwf.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/53_ffts/fftw3_fourwf.F90	2011-01-05 17:22:28.000000000 +0000
+@@ -432,8 +432,8 @@
+  end do
+ 
+  cplex=0; istwf_k=1; option=3
+-!$   call sg_fftrisc(cplex,dum_denpot,fofgin,dum_fofgin,fofr,gbound,gbound,istwf_k,dum_gvec,gvec,&
+-!$&    mgfft,ngfft,npwwfn,npwwfn,ldx,ldy,ldz,option,weight)
++!%   call sg_fftrisc(cplex,dum_denpot,fofgin,dum_fofgin,fofr,gbound,gbound,istwf_k,dum_gvec,gvec,&
++!%&    mgfft,ngfft,npwwfn,npwwfn,ldx,ldy,ldz,option,weight)
+ 
+  allocate(ftarr(2,ldx,ldy,ldz))
+ !This call gives weird results for R-->G, while G-->R is ok!!!!
+diff -Naur abinit-6.4.2.bak/src/53_ffts/sphere_fft.F90 abinit-6.4.2/src/53_ffts/sphere_fft.F90
+--- src/53_ffts/sphere_fft.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/53_ffts/sphere_fft.F90	2011-01-09 11:05:00.000000000 +0000
+@@ -99,7 +99,7 @@
+ !ENDDEBUG
+ 
+ !Insert cg into cfft with extra 0 s around outside:
+-!$OMP PARALLEL DO PRIVATE(i1,i2,i3) SHARED(cfft,ndat,n1,n2,n3)
++!%$OMP PARALLEL DO PRIVATE(i1,i2,i3) SHARED(cfft,ndat,n1,n2,n3)
+ !do i2=1,nd2proc*ndat
+ !do i3=1,n3
+ !do i1=1,n1
+@@ -108,7 +108,7 @@
+ !end do
+ !end do
+ !end do
+-!$OMP END PARALLEL DO
++!%$OMP END PARALLEL DO
+  cfft(:,:,:,:)=zero
+ !$OMP PARALLEL DO PRIVATE(i1,i2,i3,idat,ipw) SHARED(cfft,cg,kg_k,ndat,npw)
+ !write(6,*)'In sphere fft,i1,i2,i3,ipw,cfft=cg'
+@@ -229,7 +229,7 @@
+ !ENDDEBUG
+ 
+ !Insert cg into cfft with extra 0 s around outside:
+-!$OMP PARALLEL DO PRIVATE(i1,i2,i3) SHARED(cfft,ndat,n1,n2,n3)
++!%$OMP PARALLEL DO PRIVATE(i1,i2,i3) SHARED(cfft,ndat,n1,n2,n3)
+ !do i2=1,nd2proc*ndat
+ !do i3=1,n3
+ !do i1=1,n1
+@@ -238,7 +238,7 @@
+ !end do
+ !end do
+ !end do
+-!$OMP END PARALLEL DO
++!%$OMP END PARALLEL DO
+  cfft(:,:,:,:)=zero
+ !$OMP PARALLEL DO PRIVATE(i1,i2,i3,idat,ipw) SHARED(cfft,cg,kg_k,ndat,npw)
+ !write(6,*)'In sphere fft,i1,i2,i3,ipw,cfft=cg'
+diff -Naur abinit-6.4.2.bak/src/57_iovars/invars2.F90 abinit-6.4.2/src/57_iovars/invars2.F90
+--- src/57_iovars/invars2.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/57_iovars/invars2.F90	2011-01-05 16:51:21.000000000 +0000
+@@ -72,7 +72,7 @@
+  use defs_datatypes
+  use defs_abitypes
+ 
+- !$use m_gwdefs,      only: gw_sctype_from_name
++ !%use m_gwdefs,      only: gw_sctype_from_name
+ 
+ !This section has been created automatically by the script Abilint (TD).
+ !Do not modify the following lines by hand.
+diff -Naur abinit-6.4.2.bak/src/59_io_mpi/hdr_io.F90 abinit-6.4.2/src/59_io_mpi/hdr_io.F90
+--- src/59_io_mpi/hdr_io.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/59_io_mpi/hdr_io.F90	2011-01-05 08:53:52.000000000 +0000
+@@ -209,7 +209,7 @@
+  use defs_basis
+  use defs_abitypes
+  use m_errors
+-!$use m_header
++!%use m_header
+ 
+  use m_io_tools, only : flush_unit
+ 
+diff -Naur abinit-6.4.2.bak/src/64_atompaw/m_lmn_indices.F90 abinit-6.4.2/src/64_atompaw/m_lmn_indices.F90
+--- src/64_atompaw/m_lmn_indices.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/64_atompaw/m_lmn_indices.F90	2011-01-05 08:54:53.000000000 +0000
+@@ -251,10 +251,10 @@
+    indklmn(6,klmn)=jlm                       ! jlm
+ 
+    ! TODO additional mapping (At present it is in klmntomn)
+-   !$indklmn(7 ,klmn)=indlmn(2,ilmn)+il+1  ! im
+-   !$indklmn(8 ,klmn)=indlmn(2,jlmn)+jl+1  ! jm
+-   !$indklmn(9 ,klmn)=indlmn(3,ilmn)       ! in
+-   !$indklmn(10,klmn)=indlmn(3,jlmn)       ! jn
++   !%indklmn(7 ,klmn)=indlmn(2,ilmn)+il+1  ! im
++   !%indklmn(8 ,klmn)=indlmn(2,jlmn)+jl+1  ! jm
++   !%indklmn(9 ,klmn)=indlmn(3,ilmn)       ! in
++   !%indklmn(10,klmn)=indlmn(3,jlmn)       ! jn
+ 
+    if (ilm==jlm) klm_diag(klmn)=1
+   end do
+@@ -329,8 +329,8 @@
+ 
+      il= indlmn(1,ilmn); iln=indlmn(5,ilmn)
+      jl= indlmn(1,jlmn); jln=indlmn(5,jlmn)
+-     !$in = indklmn(9 ,klmn) !=indlmn(3,ilmn)
+-     !$jn = indklmn(10,klmn) !=indlmn(3,jlmn)
++     !%in = indklmn(9 ,klmn) !=indlmn(3,ilmn)
++     !%jn = indklmn(10,klmn) !=indlmn(3,jlmn)
+ 
+      in = indlmn(3,ilmn)
+      jn = indlmn(3,jlmn)
+diff -Naur abinit-6.4.2.bak/src/65_nonlocal/m_hamiltonian.F90 abinit-6.4.2/src/65_nonlocal/m_hamiltonian.F90
+--- src/65_nonlocal/m_hamiltonian.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/65_nonlocal/m_hamiltonian.F90	2011-01-05 16:38:16.000000000 +0000
+@@ -142,7 +142,7 @@
+    ! "V": all eigenvalues in the half-open interval (VL,VU] will be found.
+    ! "I": the IL-th through IU-th eigenvalues will be found.
+ 
+-  !$character(len=fnlen) :: fname
++  !%character(len=fnlen) :: fname
+   ! The name of the file storing the eigenvectors and eigenvalues (only if jobz="V")
+ 
+  end type ddiago_ctl_type
+@@ -391,7 +391,7 @@
+  gs_hamk%nfft       =nfft
+  gs_hamk%ngfft(:)   =ngfft(:)
+  gs_hamk%nloalg(:)  =nloalg(:)
+- !$gs_hamk%matblk=nloalg(4); if (nloalg(1)>0) gs_hamk%matblk=natom
++ !%gs_hamk%matblk=nloalg(4); if (nloalg(1)>0) gs_hamk%matblk=natom
+  gs_hamk%nspinor    =nspinor
+  gs_hamk%ntypat     =ntypat
+ 
+diff -Naur abinit-6.4.2.bak/src/65_psp/psp5lo.F90 abinit-6.4.2/src/65_psp/psp5lo.F90
+--- src/65_psp/psp5lo.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/65_psp/psp5lo.F90	2011-01-05 16:38:46.000000000 +0000
+@@ -111,8 +111,8 @@
+  call ctrap(mmax,work,al,result)
+ !Do integral from r(mmax) to infinity
+ !compute decay length lambda at r(mmax)
+-!$\lambda=-\log((rad(im1)*vloc(im1)+zion)$/ &
+-!$(rad(imat)*vloc(imat)+zion))/(rad(im1)-rad(imat))$
++!%\lambda=-\log((rad(im1)*vloc(im1)+zion)$/ &
++!%(rad(imat)*vloc(imat)+zion))/(rad(im1)-rad(imat))$
+ !rmtoin=$(rad(mmax)*vloc(mmax)+zion)*(rad(mmax)+1.d0/\lambda)/\lambda$
+ !Due to inability to fit exponential decay to r*V(r)+Zv
+ !in tail, NO TAIL CORRECTION IS APPLIED
+diff -Naur abinit-6.4.2.bak/src/66_paw/m_paw_slater.F90 abinit-6.4.2/src/66_paw/m_paw_slater.F90
+--- src/66_paw/m_paw_slater.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/66_paw/m_paw_slater.F90	2011-01-05 08:57:00.000000000 +0000
+@@ -352,7 +352,7 @@
+      Slatang3l(klm)%sggselect = 0
+ 
+      do ilsl=lslat_min,lslat_max
+-     !$do ilsl=lslat_min,lslat_max,2
++     !%do ilsl=lslat_min,lslat_max,2
+        lsl = ilsl-1
+        do ilc=1,lc_max
+          lc = ilc-1
+@@ -695,7 +695,7 @@
+  end if
+ 
+  indlmn  => Psps%indlmn(1:6,1:lmn_size,itypat)
+- !$indlmn => Pawtab%indklmn(1:6,1:lmn2_size) !TODO indlmn should be in Pawtab!!!
++ !%indlmn => Pawtab%indklmn(1:6,1:lmn2_size) !TODO indlmn should be in Pawtab!!!
+  indklmn => Pawtab%indklmn(1:6,1:lmn2_size)
+ 
+  allocate(kln2ln(6,ln2_size))
+@@ -874,7 +874,7 @@
+  lm2_size  = lm_size*(lm_size+1)/2
+ 
+  indlmn  => Psps%indlmn(1:6,1:lmn_size,itypat) !this is the reason why we still need itypat
+- !$indlmn => Pawtab%indlmn(1:6,1:lmn_size)     !TODO indlmn should be in Pawtab!!!
++ !%indlmn => Pawtab%indlmn(1:6,1:lmn_size)     !TODO indlmn should be in Pawtab!!!
+  indklmn => Pawtab%indklmn(1:6,1:lmn2_size)
+ 
+  ! * Setup of useful tables.
+@@ -1327,7 +1327,7 @@
+  !
+  ! Useful table for looping.
+  indlmn  => Psps%indlmn(1:6,1:lmn_size,itypat)
+- !$indlmn => Pawtab%indlmn(1:6,1:lmn2_size) !TODO indlmn should be in Pawtab!!!
++ !%indlmn => Pawtab%indlmn(1:6,1:lmn2_size) !TODO indlmn should be in Pawtab!!!
+  indklmn => Pawtab%indklmn(1:6,1:lmn2_size)
+ 
+  allocate(kln2ln(6,ln2_size))
+@@ -1822,7 +1822,7 @@
+  slat_intg=zero
+  if (Slatrad4(slt_idx)%nintgl>0) then
+    do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max
+-   !$do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max,2
++   !%do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max,2
+      isel = Slatrad4(slt_idx)%intgl_select(ilsum)
+      if (isel/=0) then
+        sltL_ijkl = Slatrad4(slt_idx)%intgl(isel)
+diff -Naur abinit-6.4.2.bak/src/66_paw/m_paw_toolbox.F90 abinit-6.4.2/src/66_paw/m_paw_toolbox.F90
+--- src/66_paw/m_paw_toolbox.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/66_paw/m_paw_toolbox.F90	2011-01-05 08:57:50.000000000 +0000
+@@ -1102,8 +1102,8 @@
+ !MG This is the coding PRESENTLY used in pawdenpot but the commented code should be the correct one
+ !MT: don't agreee for nspden (there is no dependance of vH^(1) with nspden)
+     allocate(Paw_an(iat)%vh1(cplex*Pawtab(itypat)%mesh_size,1,1))
+-    !$allocate(Paw_an(iat)%vh1 (cplex*Pawtab(itypat)%mesh_size,lm_size,nspden))
+-    !$allocate(Paw_an(iat)%vht1(cplex*Pawtab(itypat)%mesh_size,lm_size,nspden))
++    !%allocate(Paw_an(iat)%vh1 (cplex*Pawtab(itypat)%mesh_size,lm_size,nspden))
++    !%allocate(Paw_an(iat)%vht1(cplex*Pawtab(itypat)%mesh_size,lm_size,nspden))
+    end if
+   end if
+ 
+diff -Naur abinit-6.4.2.bak/src/66_paw/pawmkaewf.F90 abinit-6.4.2/src/66_paw/pawmkaewf.F90
+--- src/66_paw/pawmkaewf.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/66_paw/pawmkaewf.F90	2011-01-04 17:16:03.000000000 +0000
+@@ -99,7 +99,7 @@
+ 
+ #if defined HAVE_ETSF_IO
+  use etsf_io
+- !use etsf_io_file,     only : etsf_io_file_merge
++ use etsf_io_file,     only : etsf_io_file_merge
+  use m_abi_etsf,       only : abi_etsf_dims_init
+ #endif
+ 
+@@ -890,7 +890,7 @@
+      out_file="test_merge"
+      write(msg,'(2a)')'Master node is merging NETCDF partial files into: ',TRIM(out_file)
+      call wrtout(std_out, msg,'COLL')
+-     !$call etsf_io_file_merge(out_file,merge_files,lstat,Error)
++     call etsf_io_file_merge(out_file,merge_files,lstat,Error)
+      if (.not.lstat) goto 1000
+    end if
+    call xbarrier_mpi(spaceComm)
+diff -Naur abinit-6.4.2.bak/src/66_wfs/envlop.F90 abinit-6.4.2/src/66_wfs/envlop.F90
+--- src/66_wfs/envlop.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/66_wfs/envlop.F90	2011-01-05 16:40:07.000000000 +0000
+@@ -67,7 +67,7 @@
+ 
+ ! *************************************************************************
+ 
+-!$(k+G)^2$ evaluated using metric and kpoint
++!%$(k+G)^2$ evaluated using metric and kpoint
+  gsq(i1,i2,i3)=gmet(1,1)*(kpoint(1)+dble(i1))**2+&
+ & gmet(2,2)*(kpoint(2)+dble(i2))**2+&
+ & gmet(3,3)*(kpoint(3)+dble(i3))**2+&
+@@ -78,7 +78,7 @@
+ !cutoff function
+  fcut(gs)=(1.0_dp-(gs/cutoff))**power
+ 
+-!$(k+G)^2$ cutoff from $(1/2)(2 Pi (k+G))^2 = ecut$
++!%$(k+G)^2$ cutoff from $(1/2)(2 Pi (k+G))^2 = ecut$
+  kpgsqc=ecut/(2.0_dp*pi**2)
+  cutoff = kpgsqc
+ 
+diff -Naur abinit-6.4.2.bak/src/66_wfs/m_bands_sym.F90 abinit-6.4.2/src/66_wfs/m_bands_sym.F90
+--- src/66_wfs/m_bands_sym.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/66_wfs/m_bands_sym.F90	2011-01-05 16:43:49.000000000 +0000
+@@ -172,7 +172,7 @@
+   ! Number of states in each degenerate subspace. Cannot be larger that nclass provided
+   ! that no accidental degeneracy occurs.
+ 
+-  !$integer,pointer :: class_ids(:,:)   SET2NULL 
++  !%integer,pointer :: class_ids(:,:)   SET2NULL 
+   ! class_ids(2,nclass)
+   ! (1,icl) = index of the first symmetry of class icl
+   ! (2,icl) = index of the last symmetry of class icl
+@@ -498,12 +498,12 @@
+      MSG_ERROR("Not coded")
+ 
+      ! Read little groups from the external database.
+-     !$call init_groupk_from_file(Lgrp,spgroup,lgroup_fname,ierr)
++     !%call init_groupk_from_file(Lgrp,spgroup,lgroup_fname,ierr)
+ 
+      ! Save the irreducible representations in BSym.
+      ! Reorder symmetries such that they correspond to the Bilbao database.
+-     !$allocate(BSym%Ref_irreps(BSym%nclass)) 
+-     !$call copy_irrep(Irreps, BSym%Ref_irreps)
++     !%allocate(BSym%Ref_irreps(BSym%nclass)) 
++     !%call copy_irrep(Irreps, BSym%Ref_irreps)
+ 
+    else 
+      write(msg,'(7a)')&
+@@ -769,8 +769,8 @@
+ 
+  deallocate(sgk,tr_sgk,elements_idx)
+ 
+- !$allocate(BSym%irrep2b(0:BSym%nclass))
+- !$call nullify_coeff(BSym%irrep2b)
++ !%allocate(BSym%irrep2b(0:BSym%nclass))
++ !%call nullify_coeff(BSym%irrep2b)
+  !
+  ! 1) Allocate space for the irreducible representations.
+ 
+diff -Naur abinit-6.4.2.bak/src/66_wfs/m_wfs.F90 abinit-6.4.2/src/66_wfs/m_wfs.F90
+--- src/66_wfs/m_wfs.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/66_wfs/m_wfs.F90	2011-01-05 16:42:27.000000000 +0000
+@@ -95,7 +95,7 @@
+    ! The boundary of the basis sphere of G vectors at a given k point.
+    ! for use in improved zero padding of ffts in 3 dimensions.
+ 
+-   !$real(dp) :: kpoint(3)
++   !%real(dp) :: kpoint(3)
+ 
+    !real(dp),pointer :: ph1d(:,:)    SET2NULL
+    ! ph1d(2,3*(2*mgfft+1)*natom) 
+@@ -113,20 +113,20 @@
+    ! fnlylm(npw,dim_fnlylm,lmnmax,ntypat)
+    ! nonlocal form factors
+ 
+-   !$real(dp),pointer :: ylm(:,:)        SET2NULL
++   !%real(dp),pointer :: ylm(:,:)        SET2NULL
+    ! ylm(npw,mpsang**2*useylm)
+    ! Real spherical harmonics for each G
+ 
+-   !$integer :: ngfft(18)
++   !%integer :: ngfft(18)
+    ! Information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft.
+ 
+-   !$integer :: mgfft  
++   !%integer :: mgfft  
+    ! MAXVAL(ngfft(1:3)), used to dimension some arrays.
+ 
+-   !$integer :: nfftot 
++   !%integer :: nfftot 
+    ! PRODUCT(ngfft(1:3)), ie the total number of FFT points. 
+ 
+-   !$integer :: nfft   
++   !%integer :: nfft   
+    ! The number of points treated by this node.
+ 
+  end type kdata_t
+@@ -229,7 +229,7 @@
+   integer :: lmnmax
+   integer :: mband              ! MAX(nband)
+   integer :: mgfft              ! Maximum size of 1D FFTs
+-  !$integer :: mpsang
++  !%integer :: mpsang
+   integer :: natom
+   integer :: nfft               ! Number of FFT points treated by this processor
+   integer :: nfftot             ! Total number of points in the FFT grid
+@@ -247,7 +247,7 @@
+   integer :: my_rank            ! The rank of my processor inside the MPI communicator comm.
+   integer :: nproc              ! The number of processors in MPI comm.
+   
+-  !$ integer :: cplex
++  !% integer :: cplex
+   ! cplex= if 1 , wavefunctions are real, if 2 they are complex
+   ! In systems with time-reversal and spatial inversion, wavefunctions are real.
+   ! One might use this to reduce memory in wave_t.
+@@ -311,7 +311,7 @@
+   ! nlmn_type(ntypat)
+   ! Number of (n,l,m) channels for each type of atom. Only for PAW.
+ 
+-  !$integer,pointer :: typat(:)  SET2NULL
++  !%integer,pointer :: typat(:)  SET2NULL
+   ! typat(natom) 
+   ! The type of each atom.
+ 
+@@ -343,7 +343,7 @@
+   ! keep(mband,nkibz,nsppol)
+   ! Storage strategy: keep or not keep calculated u(r) in memory.
+ 
+-  !$type(mpicomm_t),pointer :: comm_spin(:)
++  !%type(mpicomm_t),pointer :: comm_spin(:)
+   ! comm(0:nsppol)
+   ! MPI communicators. 
+   ! comm(1), comm(2) are the MPI communicators of the nodes treating the different spins.
+@@ -978,7 +978,7 @@
+ 
+ #if 0 
+    ! TODO enable this call when k-centered G-spheres are used.
+-   !$if (wfd_ihave_ug(Wfd,0,ik_ibz,0)) then
++   !%if (wfd_ihave_ug(Wfd,0,ik_ibz,0)) then
+    call kdata_init(Wfd%Kdata(ik_ibz),Cryst,Psps,kpoint,istwf_k,ngfft,Wfd%MPI_enreg,kg_k=Wfd%gvec)
+    !endif
+ #else
+@@ -1735,12 +1735,12 @@
+  dimffnl =  Wfd%Kdata(ik_ibz)%dim_fnlylm 
+  ph3d    => Wfd%Kdata(ik_ibz)%ph3d
+  ffnl    => Wfd%Kdata(ik_ibz)%fnlylm
+- !$phkxred => Wfd%Kdata(ik_ibz)%phkxred
++ !%phkxred => Wfd%Kdata(ik_ibz)%phkxred
+ 
+ ! Compute (k+G) vectors
+  nkpg=0
+- !$if (choice==3.or.choice==2.or.choice==23) nkpg=3*Wfd%nloalg(5)
+- !$if (choice==4.or.choice==24) nkpg=9*Wfd%nloalg(5)
++ !%if (choice==3.or.choice==2.or.choice==23) nkpg=3*Wfd%nloalg(5)
++ !%if (choice==4.or.choice==24) nkpg=9*Wfd%nloalg(5)
+  allocate(kpg(npw_k,nkpg)); if (nkpg>0) call mkkpg(kg_k,kpg,kpoint,nkpg,npw_k)
+ 
+  ! Allocate and compute the arrays phkxred and ph3d
+@@ -3793,11 +3793,11 @@
+    end if
+  end do
+ 
+- !$mcg = npw_k*Wfd%nspinor*(Wfd%my_maxb-Wfd%my_minb+1)
++ !%mcg = npw_k*Wfd%nspinor*(Wfd%my_maxb-Wfd%my_minb+1)
+  ABI_CHECK(SIZE(cg,DIM=2)>=(mcg+ikg),"cg too small")
+ 
+  do band=1,Wfd%nband(ik_ibz,spin) 
+-   !$cg_bpad=npw_k*Wfd%nspinor*(band-Wfd%my_minb)
++   !%cg_bpad=npw_k*Wfd%nspinor*(band-Wfd%my_minb)
+    do ispinor=1,Wfd%nspinor
+      cg_spad=(ispinor-1)*npw_k
+      gw_spad=(ispinor-1)*Wfd%npwwfn
+diff -Naur abinit-6.4.2.bak/src/67_common/hartre1.F90 abinit-6.4.2/src/67_common/hartre1.F90
+--- src/67_common/hartre1.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/67_common/hartre1.F90	2011-01-05 17:31:18.000000000 +0000
+@@ -213,7 +213,7 @@
+ !numerator:   $1/2$ [double counting] $\times  4\pi$ [Poisson eq] $= 2\pi$
+ !denominator: $2\pi$ [reciprocal lattice vectors] squared         $= (2\pi)^2$
+ !gives the same result as real space evaluation via
+-!$\frac{1}{2}\int n(\vec r)*V_{H}(\vec r) d^3r$
++!%$\frac{1}{2}\int n(\vec r)*V_{H}(\vec r) d^3r$
+  sum(:)=sum(:)*ucvol*2*pi/(2*pi)**2
+ !MF
+ 
+diff -Naur abinit-6.4.2.bak/src/67_common/ks_ddiago.F90 abinit-6.4.2/src/67_common/ks_ddiago.F90
+--- src/67_common/ks_ddiago.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/67_common/ks_ddiago.F90	2011-01-05 17:32:44.000000000 +0000
+@@ -199,7 +199,7 @@
+  isppol  = Diago_ctl%isppol
+  kpoint  = Diago_ctl%kpoint
+  istwf_k = Diago_ctl%istwf_k
+-!$nband_k = Diago_ctl%nband_k
++!%nband_k = Diago_ctl%nband_k
+  npw_k   = Diago_ctl%npw_k
+  nloalg  = Diago_ctl%nloalg
+  ecut    = Diago_ctl%ecut
+@@ -405,7 +405,7 @@
+    end do
+  end if
+ 
+-!$call finalize_hamiltonian(gs_hamk,isppol,npw_k,istwfk,kpoint,Paw_ij)
++!%call finalize_hamiltonian(gs_hamk,isppol,npw_k,istwfk,kpoint,Paw_ij)
+ 
+ !Prepare the call to getghc.
+  ndat=1; lambda=zero; type_calc=0         ! For applying the whole Hamiltonian
+diff -Naur abinit-6.4.2.bak/src/67_common/m_coulombian.F90 abinit-6.4.2/src/67_common/m_coulombian.F90
+--- src/67_common/m_coulombian.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/67_common/m_coulombian.F90	2011-01-05 17:28:43.000000000 +0000
+@@ -1950,7 +1950,7 @@
+  real(dp) :: F3
+ !************************************************************************
+ 
+- !$F3(z)=z*\sin(qpg_para_*z)/\sqrt(rcut^2+z^2)$
++ !%$F3(z)=z*\sin(qpg_para_*z)/\sqrt(rcut^2+z^2)$
+  F3=xx*SIN(qpg_para_*xx)/SQRT(rcut_**2+xx**2)
+ 
+ end function F3
+@@ -2039,7 +2039,7 @@
+  real(dp) :: k0,rho,arg
+ !************************************************************************
+ 
+- !$K0cos(y)=K0(\rho*|qpg_z|)*COS(x.qpg_x+y*qpg_y)$
++ !%$K0cos(y)=K0(\rho*|qpg_z|)*COS(x.qpg_x+y*qpg_y)$
+  rho=SQRT(xx_**2+yy**2) ; arg=qpg_para_*rho
+  call CALCK0(arg,k0,1)
+  K0cos=k0*COS(qpgx_*xx_+qpgy_*yy)
+@@ -2065,7 +2065,7 @@
+  real(dp) :: quad
+ !************************************************************************
+ 
+- !$ K0cos_dy(x)=\int_{-b/2}^{b/2} K0(|qpg_z|\rho)cos(x.qpg_x+y.qpg_y)dy$
++ !%$K0cos_dy(x)=\int_{-b/2}^{b/2} K0(|qpg_z|\rho)cos(x.qpg_x+y.qpg_y)dy$
+  xx_=xx 
+  call quadrature(K0cos,-hb_,+hb_,qopt_,quad,ierr,ntrial_,accuracy_,npts_)
+  if (ierr/=0) then 
+diff -Naur abinit-6.4.2.bak/src/67_common/m_io_kss.F90 abinit-6.4.2/src/67_common/m_io_kss.F90
+--- src/67_common/m_io_kss.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/67_common/m_io_kss.F90	2011-01-05 17:30:48.000000000 +0000
+@@ -706,7 +706,7 @@
+ 
+  write(msg,'(3a)')ch10,' Opening file for KS structure output: ',TRIM(filekss)
+  call wrtout(std_out,msg,'COLL')
+- !$call wrtout(ab_out,msg,'COLL')
++ !%call wrtout(ab_out,msg,'COLL')
+ 
+  write(msg,'(a,i6)') ' number of Gamma centered plane waves ',npwkss
+  call wrtout(std_out,msg,'COLL') 
+diff -Naur abinit-6.4.2.bak/src/67_common/m_oscillators.F90 abinit-6.4.2/src/67_common/m_oscillators.F90
+--- src/67_common/m_oscillators.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/67_common/m_oscillators.F90	2011-01-05 17:29:29.000000000 +0000
+@@ -725,7 +725,7 @@
+  Osc%iq_bz = iq_bz
+  Osc%q_bz  = Qmesh%bz(:,iq_bz)
+ 
+- !$call get_BZ_item(Qmesh,Osc%iq_bz,Osc%q_bz,Osc%iq_ibz,Osc%isym_q,Osc%itim_q,Osc%ph_mqt)
++ !%call get_BZ_item(Qmesh,Osc%iq_bz,Osc%q_bz,Osc%iq_ibz,Osc%isym_q,Osc%itim_q,Osc%ph_mqt)
+ 
+  ! Get index of k-q = BZ(k-q) + g0. Note that k-q might fall outside the first BZ.
+  call get_BZ_diff(Kmesh,Osc%k_bz,Osc%q_bz,Osc%ikmq_bz,Osc%g0,nfound) 
+diff -Naur abinit-6.4.2.bak/src/67_common/outkss.F90 abinit-6.4.2/src/67_common/outkss.F90
+--- src/67_common/outkss.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/67_common/outkss.F90	2011-01-05 17:23:51.000000000 +0000
+@@ -132,7 +132,7 @@
+ #endif
+ 
+  use m_io_tools,         only : get_unit
+- !$use m_numeric_tools,  only : bisect
++ !%use m_numeric_tools,  only : bisect
+  use m_gsphere,          only : merge_and_sort_kg, table_gbig2kg, get_kg
+  use m_io_kss,           only : write_kss_wfgk, write_kss_header, k2gamma_centered
+  use m_hamiltonian,      only : ddiago_ctl_type, init_ddiago_ctl
+@@ -461,7 +461,7 @@
+      EXIT
+    end if
+  end do
+-!$ ishm=bisect(shlim,npwkss)
++!% ishm=bisect(shlim,npwkss)
+ 
+  if (shlim(ishm)/=npwkss) then
+    nrst1=shlim(ishm)
+diff -Naur abinit-6.4.2.bak/src/68_gw/bloch_interp.F90 abinit-6.4.2/src/68_gw/bloch_interp.F90
+--- src/68_gw/bloch_interp.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/bloch_interp.F90	2011-01-05 16:23:50.000000000 +0000
+@@ -466,7 +466,7 @@
+ 
+      call get_BZ_item(Kmesh,ik2_bz,kk2,ik2_ibz,k2_sym,k2_tim,k2_eimkt,k2_umklp)
+ 
+-     !$ik2_ibz = Kmesh%tab(ik2_bz)
++     !%ik2_ibz = Kmesh%tab(ik2_bz)
+      nband_k2 = Wfd%nband(ik2_ibz,spin)
+ 
+      inv_k2_sym = toinv(1,k2_sym)
+@@ -567,19 +567,19 @@
+    write(*,*)"with_sym",with_sym," without_sym",without_sym
+ #endif
+ 
+-   !$if (Wfd%usepaw==1) then
+-   !$  allocate(Cp_k1 (Wfd%natom,nspinor))
+-   !$  call cprj_alloc(Cp_k1,0,Wfd%nlmn_atm)
+-
+-   !$  allocate(Cp_k2 (Wfd%natom,nspinor))
+-   !$  call cprj_alloc(Cp_k2,0,Wfd%nlmn_atm)
+-   !$end if
++   !%if (Wfd%usepaw==1) then
++   !%  allocate(Cp_k1 (Wfd%natom,nspinor))
++   !%  call cprj_alloc(Cp_k1,0,Wfd%nlmn_atm)
++
++   !%  allocate(Cp_k2 (Wfd%natom,nspinor))
++   !%  call cprj_alloc(Cp_k2,0,Wfd%nlmn_atm)
++   !%end if
+ 
+    call wrtout(std_out,"Using version without symmetries","COLL")
+    do ik2_bz=1,Kmesh%nbz
+      call get_BZ_item(Kmesh,ik2_bz,kk2,ik2_ibz,k2_sym,k2_tim,k2_eimkt,k2_umklp,k2_isirred)
+ 
+-     !$ik2_ibz = Kmesh%tab(ik2_bz)
++     !%ik2_ibz = Kmesh%tab(ik2_bz)
+ 
+      nband_k2 = Wfd%nband(ik2_ibz,spin)
+ 
+@@ -597,7 +597,7 @@
+ 
+        do ik1_bz=1,ik2_bz 
+ 
+-         !$ik1_ibz = Kmesh%tab(ik1_bz)
++         !%ik1_ibz = Kmesh%tab(ik1_bz)
+          call get_BZ_item(Kmesh,ik1_bz,kk1,ik1_ibz,k1_sym,k1_tim,k1_eimkt,k1_umklp,k1_isirred)
+ 
+          nband_k1 = Wfd%nband(ik1_ibz,spin)
+@@ -616,7 +616,7 @@
+            end if
+ 
+            if (Wfd%usepaw==1) then
+-             !$call paw_symcprj(ik1_bz,nspinor,1,Cryst,Kmesh,Psps,Pawtab,Pawang,Cp_k1) 
++             !%call paw_symcprj(ik1_bz,nspinor,1,Cryst,Kmesh,Psps,Pawtab,Pawang,Cp_k1) 
+              Cp_k1 => Pku_bz(:,:,band1,ik1_bz)
+              ovlp_paw = paw_overlap(Cp_k1,Cp_k2,typat_sort,Pawtab) ! Be careful as Cp are always sorted.
+              blk_ovlp = blk_ovlp + DCMPLX(ovlp_paw(1),ovlp_paw(2))
+@@ -671,8 +671,8 @@
+        end do
+      end do
+      deallocate(Pku_bz)
+-     !$call cprj_free(Cp_k1); deallocate(Cp_k1)
+-     !$call cprj_free(Cp_k2); deallocate(Cp_k2)
++     !%call cprj_free(Cp_k1); deallocate(Cp_k1)
++     !%call cprj_free(Cp_k2); deallocate(Cp_k2)
+    end if
+    !
+    ! ========================================
+@@ -718,8 +718,8 @@
+    !
+ 
+    !istwfk=1; nsppol_=1; nkibz_=1: nband_(:,:)=blc_size
+-   !$call wfd_init(Wlcd,Cryst,Pawtab,Psps,keep_ur,Wfd%paral_kgb,Wfd%npwwfn,blc_size,nband_,nkibz_,nsppol_,bks_mask,&
+-   !$&  Wfd%nspden,Wfd%nspinor,istwfk,kibz,ngfft,mg0,gvec,nloalg,comm)
++   !%call wfd_init(Wlcd,Cryst,Pawtab,Psps,keep_ur,Wfd%paral_kgb,Wfd%npwwfn,blc_size,nband_,nkibz_,nsppol_,bks_mask,&
++   !%&  Wfd%nspden,Wfd%nspinor,istwfk,kibz,ngfft,mg0,gvec,nloalg,comm)
+ 
+    allocate(blc_ug(Wfd%npwwfn*nspinor,blc_size))
+    allocate(sum_ur(Wfd%nfft*nspinor))
+@@ -799,7 +799,7 @@
+ 
+      write(blc_unt,rec=1)reclen,Wfd%npwwfn,nspinor,blc_size,spin
+      do blc=1,blc_size 
+-       !$write(blc_unt,rec=1+blc)blc_ug(:,blc)
++       !%write(blc_unt,rec=1+blc)blc_ug(:,blc)
+      end do
+ 
+      close(unit=blc_unt)
+@@ -1274,11 +1274,11 @@
+ 
+  ! Recalculate new occupation numbers and new Fermi level.
+  ! FIXME Have to calculate new occ for semiconductors
+- !$call update_occ(ioBSt,Dtset%fixmom,Dtset%stmbias,Dtset%prtvol)
++ !%call update_occ(ioBSt,Dtset%fixmom,Dtset%stmbias,Dtset%prtvol)
+ 
+  call bst_write_bands(ioBSt,Cryst%gmet,"interpolated",ierr)
+ 
+- !$call bst_print_fs(ioBst,Cryst,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk,"Fermi_surface",ierr)            
++ !%call bst_print_fs(ioBst,Cryst,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk,"Fermi_surface",ierr)            
+  !
+  ! ===================================
+  ! ==== Optionally Write WFK file ====
+@@ -1380,8 +1380,8 @@
+        end do
+        deallocate(umat_k)
+        
+-       !$call XGEMM('N','N',nsize,sizegw,rangeb,cone_gw,wf_ks(:,lowerb:upperb),nsize,&
+-       !$&          umat_k,rangeb,czero_gw,wf_qp(spad+1:spad+nsize,b1gw:b2gw),nsize)
++       !%call XGEMM('N','N',nsize,sizegw,rangeb,cone_gw,wf_ks(:,lowerb:upperb),nsize,&
++       !%&          umat_k,rangeb,czero_gw,wf_qp(spad+1:spad+nsize,b1gw:b2gw),nsize)
+ 
+ #if 0
+        ! * Perform conversion of the basis set.
+diff -Naur abinit-6.4.2.bak/src/68_gw/calc_density.F90 abinit-6.4.2/src/68_gw/calc_density.F90
+--- src/68_gw/calc_density.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/calc_density.F90	2011-01-05 16:15:05.000000000 +0000
+@@ -336,8 +336,8 @@
+  end if
+  call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
+ 
+-!$write(msg,'(a,f9.4)')' Renormalizing smooth charge density using nratio = ',nratio
+-!$ rhor(:,:)=nratio*rhor(:,:)
++!%write(msg,'(a,f9.4)')' Renormalizing smooth charge density using nratio = ',nratio
++!% rhor(:,:)=nratio*rhor(:,:)
+ 
+  write(msg,'(a,f9.6)')' average of density, n = ',rhoav
+  call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
+diff -Naur abinit-6.4.2.bak/src/68_gw/calc_sigc_cd.F90 abinit-6.4.2/src/68_gw/calc_sigc_cd.F90
+--- src/68_gw/calc_sigc_cd.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/calc_sigc_cd.F90	2011-01-05 16:28:20.000000000 +0000
+@@ -137,7 +137,7 @@
+      call spline(DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),nomegaer,zero,zero,rtmp_r)
+      call spline(DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),nomegaer,zero,zero,rtmp_i)
+ 
+-     !$call spline_complex( DBLE(omega(1:nomegaer)), epsrho(ig,1:nomegaer), nomegaer, zero, zero, rtmp )
++     !%call spline_complex( DBLE(omega(1:nomegaer)), epsrho(ig,1:nomegaer), nomegaer, zero, zero, rtmp )
+ 
+      do ios=1,nomega
+ 
+@@ -150,7 +150,7 @@
+        call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),rtmp_i,1,tmp_x,tmp_y)
+        rt_imag = tmp_y(1)
+ 
+-       !$call splint_complex(nomegaer,DBLE(omega(1:nomegaer)),epsrho(ig,1:nomegaer),rtmp,1,tmp_x,yfit)
++       !%call splint_complex(nomegaer,DBLE(omega(1:nomegaer)),epsrho(ig,1:nomegaer),rtmp,1,tmp_x,yfit)
+ 
+        ct=DCMPLX(rt_real,rt_imag)
+ 
+diff -Naur abinit-6.4.2.bak/src/68_gw/calc_sigc_me.F90 abinit-6.4.2/src/68_gw/calc_sigc_me.F90
+--- src/68_gw/calc_sigc_me.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/calc_sigc_me.F90	2011-01-05 16:35:54.000000000 +0000
+@@ -266,7 +266,7 @@
+  ! Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
+  jk_bz=Sigp%kptgw2bz(ikcalc)
+  call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
+- !$call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
++ !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
+  spinrot_kgw=Cryst%spinrot(:,isym_kgw)
+  !
+  ib1=minbnd
+@@ -450,10 +450,10 @@
+    !   pass it to the setup_ppmodel
+    ! * It would be possible to calculate rho(G) using Paw_pwff though. Maybe faster but
+    !   results would depend on the expression used for the matrix elements. This is safer.
+-!$   allocate(ks_rhor_paw(rho_nfftot,Dtset%nspden))
+-!$   call denfgr(Cryst%natom,Dtset%nspden,ks_nhat,Cryst%ntypat,Pawfgr,Pawfgrtab,Pawrad,KS_pawrhoij,Pawtab,&
+-!$   &   ks_rhor,ks_rhor_paw,Psps,Cryst%typat)
+-!$   deallocate(ks_rhor_paw)
++!%   allocate(ks_rhor_paw(rho_nfftot,Dtset%nspden))
++!%   call denfgr(Cryst%natom,Dtset%nspden,ks_nhat,Cryst%ntypat,Pawfgr,Pawfgrtab,Pawrad,KS_pawrhoij,Pawtab,&
++!%   &   ks_rhor,ks_rhor_paw,Psps,Cryst%typat)
++!%   deallocate(ks_rhor_paw)
+  end if ! usepaw==1
+  !
+  !
+@@ -640,8 +640,8 @@
+ #ifdef HAVE_CLIB
+    call clib_progress_bar(ik_bz,Kmesh%nbz)
+ #else
+-   !$write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
+-   !$call wrtout(std_out,msg,'PERS')
++   !%write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
++   !%call wrtout(std_out,msg,'PERS')
+ #endif
+    !
+    ! === Find the corresponding irred q-point ===
+diff -Naur abinit-6.4.2.bak/src/68_gw/calc_sigx_me.F90 abinit-6.4.2/src/68_gw/calc_sigx_me.F90
+--- src/68_gw/calc_sigx_me.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/calc_sigx_me.F90	2011-01-05 10:06:07.000000000 +0000
+@@ -237,7 +237,7 @@
+  ! === Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
+  jk_bz=Sigp%kptgw2bz(ikcalc)
+  call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
+- !$call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
++ !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
+  spinrot_kgw(:)=Cryst%spinrot(:,isym_kgw)
+  !
+  write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,&
+@@ -486,8 +486,8 @@
+ #ifdef HAVE_CLIB
+      call clib_progress_bar(ik_bz,Kmesh%nbz)
+ #else
+-     !$write(msg,'(2(a,i4),a,i3)')' calc_sigx_me : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
+-     !$call wrtout(std_out,msg,'PERS')
++     !%write(msg,'(2(a,i4),a,i3)')' calc_sigx_me : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
++     !%call wrtout(std_out,msg,'PERS')
+ #endif
+      !
+      ! * Find the corresponding irreducible q-point.
+diff -Naur abinit-6.4.2.bak/src/68_gw/cchi0.F90 abinit-6.4.2/src/68_gw/cchi0.F90
+--- src/68_gw/cchi0.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/cchi0.F90	2011-01-05 10:06:42.000000000 +0000
+@@ -417,7 +417,7 @@
+      call clib_progress_bar(ik_bz,Kmesh%nbz)
+ #else
+      write(msg,'(2(a,i4),a,i2,a,i3)')' ik = ',ik_bz,' / ',Kmesh%nbz,' is = ',is,' done by processor ',Wfd%my_rank
+-     !$call wrtout(std_out,msg,'PERS')
++     !%call wrtout(std_out,msg,'PERS')
+ #endif
+      !
+      ! * Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz.
+diff -Naur abinit-6.4.2.bak/src/68_gw/classify_bands.F90 abinit-6.4.2/src/68_gw/classify_bands.F90
+--- src/68_gw/classify_bands.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/classify_bands.F90	2011-01-05 16:28:45.000000000 +0000
+@@ -253,7 +253,7 @@
+      !tmp_sym(:,:,isym) = Cryst%symrec(:,:,isym)
+      !tmp_sym(:,:,isym) = TRANSPOSE(Cryst%symrec(:,:,isym))
+    end do
+-   !$call setsymrhoij(Cryst%rprimd,lmax,Cryst%nsym,3,Cryst%gprimd,tmp_sym,zarot)
++   !%call setsymrhoij(Cryst%rprimd,lmax,Cryst%nsym,3,Cryst%gprimd,tmp_sym,zarot)
+    !call setsymrhoij(Cryst%gprimd,lmax,Cryst%nsym,1,Cryst%rprimd,tmp_sym,zarot)
+    deallocate(tmp_sym)
+    zarot = Pawang%zarot
+diff -Naur abinit-6.4.2.bak/src/68_gw/cohsex_me.F90 abinit-6.4.2/src/68_gw/cohsex_me.F90
+--- src/68_gw/cohsex_me.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/cohsex_me.F90	2011-01-05 10:14:27.000000000 +0000
+@@ -228,7 +228,7 @@
+  ! Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
+  jk_bz=Sigp%kptgw2bz(ikcalc)
+  call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
+- !$call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
++ !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
+  spinrot_kgw=Cryst%spinrot(:,isym_kgw)
+  !
+  ib1=minbnd
+@@ -559,8 +559,8 @@
+ #ifdef HAVE_CLIB
+      call clib_progress_bar(ik_bz,Kmesh%nbz)
+ #else
+-     !$write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
+-     !$call wrtout(std_out,msg,'PERS')
++     !%write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
++     !%call wrtout(std_out,msg,'PERS')
+ #endif
+      !
+      ! === Find the corresponding irred q-point ===
+diff -Naur abinit-6.4.2.bak/src/68_gw/debug_tools.F90 abinit-6.4.2/src/68_gw/debug_tools.F90
+--- src/68_gw/debug_tools.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/debug_tools.F90	2011-01-05 16:31:28.000000000 +0000
+@@ -624,9 +624,9 @@
+    ! FIXME should check the expression in case of non zero tnons.
+    ! /***********************************************************************/
+ 
+-   !$all init_oscillator(Osc,isppol,jkbz,Kmesh,iq_bz,Qmesh,Ep%npwe,nspinor,(/ib,ib/),(/ib1,ib2/) )
+-   !$call calc_pw_oscillator(Wf_info,Cryst,Osc,MPI_enreg)
+-   !$call destroy_oscillator(Osc)
++   !%all init_oscillator(Osc,isppol,jkbz,Kmesh,iq_bz,Qmesh,Ep%npwe,nspinor,(/ib,ib/),(/ib1,ib2/) )
++   !%call calc_pw_oscillator(Wf_info,Cryst,Osc,MPI_enreg)
++   !%call destroy_oscillator(Osc)
+ 
+    do ib1=1,Ep%nbnds ! Loop over "conduction" states.
+     do ib2=1,Ep%nbnds ! Loop over "valence" states.
+diff -Naur abinit-6.4.2.bak/src/68_gw/gw_tools.F90 abinit-6.4.2/src/68_gw/gw_tools.F90
+--- src/68_gw/gw_tools.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/gw_tools.F90	2011-01-05 10:04:31.000000000 +0000
+@@ -760,7 +760,7 @@
+        if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then 
+          bbp_mask(ib1,ib2)=.TRUE.
+          if (deltaeGW_b1kmq_b2k<zero) bbp_mask(ib1,ib2)=.FALSE.  ! Only positive frequencies are needed for the Hilbert transform.
+-         !$if (use_tr .and. ib1<ib2) bbp_mask(ib1,ib2)=.FALSE. ! GAIN a factor ~2 thanks to time-reversal.
++         !%if (use_tr .and. ib1<ib2) bbp_mask(ib1,ib2)=.FALSE. ! GAIN a factor ~2 thanks to time-reversal.
+        end if
+ 
+      CASE  DEFAULT
+diff -Naur abinit-6.4.2.bak/src/68_gw/m_dyson_solver.F90 abinit-6.4.2/src/68_gw/m_dyson_solver.F90
+--- src/68_gw/m_dyson_solver.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/m_dyson_solver.F90	2011-01-05 10:04:58.000000000 +0000
+@@ -172,7 +172,7 @@
+  ! Index of this k-point in the IBZ array.
+  ikbz_gw=Sigp%kptgw2bz(ikcalc) 
+  call get_BZ_item(Kmesh,ikbz_gw,kbz_gw,jkibz,isym,itim,phase)
+- !$call get_IBZ_item(Kmesh,jkibz,kibz,wtk)
++ !%call get_IBZ_item(Kmesh,jkibz,kibz,wtk)
+ 
+  sigc=czero; dsigc=czero
+ 
+diff -Naur abinit-6.4.2.bak/src/68_gw/m_fft_prof.F90 abinit-6.4.2/src/68_gw/m_fft_prof.F90
+--- src/68_gw/m_fft_prof.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/m_fft_prof.F90	2011-01-05 10:05:26.000000000 +0000
+@@ -64,11 +64,11 @@
+  type,public :: FFT_prof_t
+    !private
+    integer :: ncalls
+-   !$integer :: ndat
++   !%integer :: ndat
+    character(len=TNAME_LEN) :: test_name   
+    real(dp) :: cpu_time
+    real(dp) :: wall_time
+-   !$integer :: ngfft(18)=-1
++   !%integer :: ngfft(18)=-1
+    real(dp),pointer :: results(:,:,:)  SET2NULL
+  end type FFT_prof_t
+ 
+diff -Naur abinit-6.4.2.bak/src/68_gw/m_gwannier.F90 abinit-6.4.2/src/68_gw/m_gwannier.F90
+--- src/68_gw/m_gwannier.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/m_gwannier.F90	2011-01-05 16:34:34.000000000 +0000
+@@ -365,8 +365,8 @@
+  call destroy_BZ_mesh_type(Kpath)
+  call bstruct_clean(KS_BSt) 
+  call bstruct_clean(QP_BSt) 
+- !$call bstruct_clean(QP_intp) 
+- !$call bstruct_clean(KS_intp) 
++ !%call bstruct_clean(QP_intp) 
++ !%call bstruct_clean(KS_intp) 
+  call hdr_clean(Hdr)
+ 
+  STOP 'myGWannier OK'
+diff -Naur abinit-6.4.2.bak/src/68_gw/m_io_screening.F90 abinit-6.4.2/src/68_gw/m_io_screening.F90
+--- src/68_gw/m_io_screening.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/m_io_screening.F90	2011-01-05 16:19:49.000000000 +0000
+@@ -66,69 +66,69 @@
+ !!
+ !! SOURCE
+ 
+-!$ type,public :: ScrHdr_type
+-!$
+-!$  !Other variables that can be added are, for the moment, commented out. 
+-!$  !Most of them are related to the Abinit implementation  and should not compare in the ETSF specs
+-!$
+-!$  !Index of the qlwl section?
+-!$  !gwcomp, gwencomp  ! Info on the extrapolar algorithm
+-!$ 
+-!$  integer :: ID           ! Matrix identifier: O if not yet defined, 1 for chi0, 
+-!$                            ! 2 for chi, 3 for epsilon, 4 for espilon^{-1}
+-!$  integer :: ikxc         ! Kxc kernel used, 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for frequency-dependent TDDFT 
+-!$  integer :: inclvkb      ! q-->0 treatment, 0 for None, 1-2 for transversal gauge, 3 for longitudinal
+-!$  integer :: headform     ! format of the SCR header
+-!$  integer :: fform        ! File format:
+-!$  integer :: gwcalctyp    ! Calculation type (G0W0, G0W, GW ...)
+-!$  integer :: nI,nJ        ! Number of spin components (rows,columns) in chi|eps^-1. (1,1) if collinear. 
+-!$                          !  The internal representation of the matrix is eps(nI*npwe,nJ*npwe) 
+-!$  integer :: nqibz        ! Number of q-points in the IBZ.
+-!$  integer :: nqlwl        ! Number of points for the treatment of the long wavelength limit.
+-!$  integer :: nomega       ! Total number of frequencies.
+-!$  integer :: nbnds_used   ! Number of bands used during the screening calculation (only for info)
+-!$  integer :: npwe         ! Number of G vectors reported on the file.
+-!$  integer :: npwwfn_used  ! Number of G vectors for wavefunctions used during the screening calculation (only for info)
+-!$  integer :: spmeth       ! Method used to approximate the delta function in the expression for Im Chi_0
+-!$  integer :: test_type    ! 0 for None, 1 for TEST-PARTICLE, 2 for TEST-ELECTRON (only for TDDFT)
+-!$  integer :: tordering    ! 0 if not defined, 1 for Time-Ordered, 2 for Advanced, 3 for Retarded.
+-!$
+-!$  real(dp) :: soenergy    ! Scissor Energy, zero if not used
+-!$  real(dp) :: spsmear     ! Smearing of the delta in case of spmeth==2
+-!$  real(dp) :: zcut        ! Imaginary shift to avoid the poles along the real axis.
+-!$
+-!$  type(Hdr_type) :: Hdr   ! The abinit header.
+-!$
+-!$!arrays
+-!$  character(len=80) :: title(2)
+-!$  ! Title describing the content of the file.
+-!$
+-!$  integer,pointer  :: gvec(:,:)                 
+-!$  ! gvec(3,npwe) 
+-!$  ! G vectors in r.l.u.
+-!$
+-!$  real(dp),pointer :: qibz(:,:)
+-!$  ! qibz(3,nqibz)
+-!$  ! q-points in r.l.u.
+-!$
+-!$  real(dp),pointer :: qlwl(:,:)                 
+-!$  ! qlwl(3,nqlwl)
+-!$  ! q-points for the long wave-length limit treatment (r.l.u)
+-!$
+-!$  complex(dpc),pointer :: lwing(:,:,:)         
+-!$  ! lwing(npwe,nomega,nqlwl)
+-!$  ! Lower wings for the different q"s -->0 
+-!$
+-!$  complex(dpc),pointer :: omega(:)             
+-!$  ! omega(nomega) 
+-!$  ! All frequencies calculated both along the real and the imaginary axis.
+-!$
+-!$  complex(dpc),pointer :: uwing(:,:,:)   
+-!$  ! uwing(npwe,nomega,nqlwl)
+-!$  ! Upper wings for the different q"s -->0 
+-!$
+-!$ end type ScrHdr_type
+-!$!!***
++!% type,public :: ScrHdr_type
++!%
++!%  !Other variables that can be added are, for the moment, commented out. 
++!%  !Most of them are related to the Abinit implementation  and should not compare in the ETSF specs
++!%
++!%  !Index of the qlwl section?
++!%  !gwcomp, gwencomp  ! Info on the extrapolar algorithm
++!% 
++!%  integer :: ID           ! Matrix identifier: O if not yet defined, 1 for chi0, 
++!%                            ! 2 for chi, 3 for epsilon, 4 for espilon^{-1}
++!%  integer :: ikxc         ! Kxc kernel used, 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for frequency-dependent TDDFT 
++!%  integer :: inclvkb      ! q-->0 treatment, 0 for None, 1-2 for transversal gauge, 3 for longitudinal
++!%  integer :: headform     ! format of the SCR header
++!%  integer :: fform        ! File format:
++!%  integer :: gwcalctyp    ! Calculation type (G0W0, G0W, GW ...)
++!%  integer :: nI,nJ        ! Number of spin components (rows,columns) in chi|eps^-1. (1,1) if collinear. 
++!%                          !  The internal representation of the matrix is eps(nI*npwe,nJ*npwe) 
++!%  integer :: nqibz        ! Number of q-points in the IBZ.
++!%  integer :: nqlwl        ! Number of points for the treatment of the long wavelength limit.
++!%  integer :: nomega       ! Total number of frequencies.
++!%  integer :: nbnds_used   ! Number of bands used during the screening calculation (only for info)
++!%  integer :: npwe         ! Number of G vectors reported on the file.
++!%  integer :: npwwfn_used  ! Number of G vectors for wavefunctions used during the screening calculation (only for info)
++!%  integer :: spmeth       ! Method used to approximate the delta function in the expression for Im Chi_0
++!%  integer :: test_type    ! 0 for None, 1 for TEST-PARTICLE, 2 for TEST-ELECTRON (only for TDDFT)
++!%  integer :: tordering    ! 0 if not defined, 1 for Time-Ordered, 2 for Advanced, 3 for Retarded.
++!%
++!%  real(dp) :: soenergy    ! Scissor Energy, zero if not used
++!%  real(dp) :: spsmear     ! Smearing of the delta in case of spmeth==2
++!%  real(dp) :: zcut        ! Imaginary shift to avoid the poles along the real axis.
++!%
++!%  type(Hdr_type) :: Hdr   ! The abinit header.
++!%
++!%!arrays
++!%  character(len=80) :: title(2)
++!%  ! Title describing the content of the file.
++!%
++!%  integer,pointer  :: gvec(:,:)                 
++!%  ! gvec(3,npwe) 
++!%  ! G vectors in r.l.u.
++!%
++!%  real(dp),pointer :: qibz(:,:)
++!%  ! qibz(3,nqibz)
++!%  ! q-points in r.l.u.
++!%
++!%  real(dp),pointer :: qlwl(:,:)                 
++!%  ! qlwl(3,nqlwl)
++!%  ! q-points for the long wave-length limit treatment (r.l.u)
++!%
++!%  complex(dpc),pointer :: lwing(:,:,:)         
++!%  ! lwing(npwe,nomega,nqlwl)
++!%  ! Lower wings for the different q"s -->0 
++!%
++!%  complex(dpc),pointer :: omega(:)             
++!%  ! omega(nomega) 
++!%  ! All frequencies calculated both along the real and the imaginary axis.
++!%
++!%  complex(dpc),pointer :: uwing(:,:,:)   
++!%  ! uwing(npwe,nomega,nqlwl)
++!%  ! Upper wings for the different q"s -->0 
++!%
++!% end type ScrHdr_type
++!%!!***
+ 
+  public :: scr_hdr_io        ! I/O of the header (read/write/echo).
+  public :: print_ScrHdr      ! Print the SCR related part of the header.
+@@ -748,7 +748,7 @@
+  Hscr%ikxc       =ikxc         
+  Hscr%inclvkb    =Ep%inclvkb      
+  Hscr%fform      =1002
+- !$Hscr%fform      =1102
++ !%Hscr%fform      =1102
+  Hscr%headform   =HSCR_LATEST_HEADFORM
+  Hscr%gwcalctyp  =Ep%gwcalctyp
+  Hscr%nI         =Ep%nI
+diff -Naur abinit-6.4.2.bak/src/68_gw/m_qparticles.F90 abinit-6.4.2/src/68_gw/m_qparticles.F90
+--- src/68_gw/m_qparticles.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/m_qparticles.F90	2011-01-05 16:22:01.000000000 +0000
+@@ -333,7 +333,7 @@
+  ! =====================
+  
+  ! Dump info on the crystal structure.
+- !$call abi_crystal_put(Cryst,filetsf)
++ !%call abi_crystal_put(Cryst,filetsf)
+ 
+  call etsf_io_low_set_write_mode(ncid,lstat,Error_data=Error)
+  if (.not.lstat) goto 1000
+@@ -500,7 +500,7 @@
+  ABI_CHECK(dimrho==0.or.dimrho==1,'dimrho must be 0 or 1')
+ 
+  ! This does not work in parallel !!?
+- !$my_rank = xcomm_rank(MPI_enreg%spaceComm)
++ !%my_rank = xcomm_rank(MPI_enreg%spaceComm)
+  call xme_init(MPI_enreg,my_rank)
+  master=0
+  !
+@@ -691,7 +691,7 @@
+ 
+        call rhoij_io(pawrhoij,unqps,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,&
+ &                    HDR_LATEST_HEADFORM,"Read",form="formatted")
+-       !$call rhoij_io(pawrhoij,std_out,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,HDR_LATEST_HEADFORM,"Echo")
++       !%call rhoij_io(pawrhoij,std_out,BSt%nsppol,BSt%nspinor,nspden,nlmn_type,Cryst%typat,HDR_LATEST_HEADFORM,"Echo")
+ 
+        deallocate(nlmn_type,typatR)
+      end if ! usepaw
+@@ -739,9 +739,9 @@
+   deallocate(full_rmtmp)
+ 
+   ! Read energies.
+-  !$call etsf_io_low_read_var(ncid,'m_lda_to_qp',m_lda_to_qp,lstat,Error_data=Error)
+-  !$if (.not.lstat) goto 1000
+-  !$BSt%eig(1:nbsc,ik,isppol)=en_tmp(1:nbsc)
++  !%call etsf_io_low_read_var(ncid,'m_lda_to_qp',m_lda_to_qp,lstat,Error_data=Error)
++  !%if (.not.lstat) goto 1000
++  !%BSt%eig(1:nbsc,ik,isppol)=en_tmp(1:nbsc)
+ 
+   if ( ALL(ngfftf(1:3)==(/n1,n2,n3/)) ) then
+    call etsf_io_low_read_var(ncid,'qp_density',rhor_out,lstat,Error_data=Error)
+diff -Naur abinit-6.4.2.bak/src/68_gw/m_screening.F90 abinit-6.4.2/src/68_gw/m_screening.F90
+--- src/68_gw/m_screening.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/m_screening.F90	2011-01-05 16:30:36.000000000 +0000
+@@ -360,7 +360,7 @@
+    end do
+ 
+    rdwr=4
+-   !$call hdr_io_int(Er%fform,Er%Hscr%Hdr,rdwr,unt)
++   !%call hdr_io_int(Er%fform,Er%Hscr%Hdr,rdwr,unt)
+  end if ! verbose>0
+ 
+ end subroutine print_epsilonm1_results 
+@@ -712,7 +712,7 @@
+ 
+  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
+ 
+- !$ if (Er%ID/=0) call reset_Epsilonm1(Er)
++ !% if (Er%ID/=0) call reset_Epsilonm1(Er)
+  Er%ID=id_required
+ 
+  if (Er%ID==Er%Hscr%ID) then
+@@ -1790,7 +1790,7 @@
+    end if
+ 
+    if (iqibz==1) then 
+-     !$vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl)  ! Use Coulomb term for q-->0
++     !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl)  ! Use Coulomb term for q-->0
+      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! TODO add treatment of non-Analytic behavior
+    else 
+      vc_sqrt => Vcp%vc_sqrt(:,iqibz)  
+@@ -2169,7 +2169,7 @@
+    !
+    ! Change the body but do not add the corrections due to the head and the wings.
+    ! since they can be obtained on the fly from eps_body and the wings of eps^{-1}.
+-   !$chi0(2:,2:,iomega) = eps_body
++   !%chi0(2:,2:,iomega) = eps_body
+  end do !iomega
+ 
+  deallocate(modg_inv,cvec)
+diff -Naur abinit-6.4.2.bak/src/68_gw/m_sigma_results.F90 abinit-6.4.2/src/68_gw/m_sigma_results.F90
+--- src/68_gw/m_sigma_results.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/m_sigma_results.F90	2011-01-05 16:32:15.000000000 +0000
+@@ -1220,7 +1220,7 @@
+  allocate(gw_corrections(2,KS_BSt%mband,KS_BSt%nkpt,KS_BSt%nsppol))
+  gw_corrections=zero
+  gw_corrections(1,:,:,:) = QP_BSt%eig - KS_BSt%eig
+- !$gw_corrections = c2r(Sr%degw)
++ !%gw_corrections = c2r(Sr%degw)
+  GWdata%gw_corrections%data4D => gw_corrections
+ 
+  call etsf_io_gwdata_put(ncid,GWdata,lstat,Error_data)
+diff -Naur abinit-6.4.2.bak/src/68_gw/print_psps.F90 abinit-6.4.2/src/68_gw/print_psps.F90
+--- src/68_gw/print_psps.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/print_psps.F90	2011-01-05 16:29:06.000000000 +0000
+@@ -182,7 +182,7 @@
+   !  If mpspso is 2, lmnmax takes into account the spin-orbit projectors,
+   !  so, it is equal to the max of lmnprojso or lnprojso, see pspheader_type
+ 
+-!$integer, pointer :: indlmn(:,:,:)
++!%integer, pointer :: indlmn(:,:,:)
+ ! indlmn(6,lmnmax,ntypat)
+ ! For each type of psp,
+ ! array giving l,m,n,lm,ln,spin for i=ln  (if useylm=0)
+diff -Naur abinit-6.4.2.bak/src/68_gw/setup_screening.F90 abinit-6.4.2/src/68_gw/setup_screening.F90
+--- src/68_gw/setup_screening.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/setup_screening.F90	2011-01-05 16:27:08.000000000 +0000
+@@ -238,10 +238,10 @@
+  call print_BZ_mesh(Kmesh,"K-mesh for the wavefunctions",std_out,Dtset%prtvol,"COLL")
+  call print_BZ_mesh(Kmesh,"K-mesh for the wavefunctions",ab_out, 0,           "COLL")
+ 
+- !$call nullify_BZ_mesh(Kmesh4test)
+- !$call make_mesh(Kmesh4test,Cryst,Dtset%kptopt,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk)
+- !$call print_BZ_mesh(Kmesh4test,"Kmesh4test",std_out,prtvol=10)
+- !$call destroy_BZ_mesh_type(Kmesh4test)
++ !%call nullify_BZ_mesh(Kmesh4test)
++ !%call make_mesh(Kmesh4test,Cryst,Dtset%kptopt,Dtset%kptrlatt,Dtset%nshiftk,Dtset%shiftk)
++ !%call print_BZ_mesh(Kmesh4test,"Kmesh4test",std_out,prtvol=10)
++ !%call destroy_BZ_mesh_type(Kmesh4test)
+  !
+  ! === Find Q-mesh, and do setup for long wavelength limit ===
+  ! * Stop if a nonzero umklapp is needed to reconstruct the BZ. In this case, indeed, 
+@@ -465,7 +465,7 @@
+ 
+  ! To write the SCR header correctly, with heads and wings, we have 
+  ! to make sure that q==0, if present, is the first q-point in the list.
+- !$has_q0=(ANY(normv(Ep%qcalc(:,:),gmet,'G')<GW_TOLQ0)) !commented to avoid problems with sunstudio12
++ !%has_q0=(ANY(normv(Ep%qcalc(:,:),gmet,'G')<GW_TOLQ0)) !commented to avoid problems with sunstudio12
+  has_q0=.FALSE.
+  do iq=1,Ep%nqcalc
+    if (normv(Ep%qcalc(:,iq),gmet,'G')<GW_TOLQ0) then
+@@ -522,7 +522,7 @@
+  !write(*,*)MAXVAL(ABS(occfact(:)-KS_BSt%occ(:))) 
+ 
+  !TODO call update_occ here
+- !$call update_occ(KS_BSt,fixmom,stmbias,Dtset%prtvol)
++ !%call update_occ(KS_BSt,fixmom,stmbias,Dtset%prtvol)
+ 
+  deallocate(doccde,eigen,npwarr)
+ 
+diff -Naur abinit-6.4.2.bak/src/68_gw/setup_sigma.F90 abinit-6.4.2/src/68_gw/setup_sigma.F90
+--- src/68_gw/setup_sigma.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/setup_sigma.F90	2011-01-05 16:33:32.000000000 +0000
+@@ -421,7 +421,7 @@
+ &  Kmesh%nibz,npwarr,Dtset%nsppol,Dtset%nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,occfact,Kmesh%wt) 
+ 
+  !TODO call update_occ here
+- !$call update_occ(KS_BSt,fixmom,stmbias,Dtset%prtvol)
++ !%call update_occ(KS_BSt,fixmom,stmbias,Dtset%prtvol)
+ 
+  ! this fails simply because in case of NSCF occ are zero
+  ! TODO outkss should calculate occ factors in case of NSCF run. 
+diff -Naur abinit-6.4.2.bak/src/68_gw/trashme.F90 abinit-6.4.2/src/68_gw/trashme.F90
+--- src/68_gw/trashme.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_gw/trashme.F90	2011-01-05 16:29:33.000000000 +0000
+@@ -911,7 +911,7 @@
+  deallocate(umat_k)
+  deallocate(wf_ks,wf_qp)
+ 
+- !$call wfd_reset_ur(Wf_info)
++ !%call wfd_reset_ur(Wf_info)
+ 
+  DBG_EXIT("COLL")
+ 
+diff -Naur abinit-6.4.2.bak/src/68_rsprc/ladielmt.F90 abinit-6.4.2/src/68_rsprc/ladielmt.F90
+--- src/68_rsprc/ladielmt.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_rsprc/ladielmt.F90	2011-01-05 09:31:41.000000000 +0000
+@@ -246,9 +246,9 @@
+ 
+ 
+ !***********************************************************************************
+-!Getting the localy averaged non-local potential                                ***
+-!$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
+-!**********************************************************************************
++! Getting the localy averaged non-local potential                                ***
++!%$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
++!***********************************************************************************
+ 
+  qphon=zero
+ 
+@@ -671,36 +671,36 @@
+ 
+ 
+ !    lavnlr(:,ispden) = -min
+-!    !$        write(filapp,*) 'orig-lavnl-0.00'
+-!    !$        call out1line(filapp,lavnlr,dtset,0,0)
+-!    !$        write(filapp,*) 'orig-lavnl-0.25'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
+-!    !$        write(filapp,*) 'orig-lavnl-0.50'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
+-!    !$        lavnlr=lavnlr*(rhor+0.0001_dp)
+-!    !$ !DEBUG
+-!    !$write(0,*) "DEBUG ================ lavnl ==========================="
+-!    !$write(0,*) "DEBUG lavnl: density analysis"
+-!    !$ !ENDDEBUG
+-!    !$        do ifft=1,dtset%nfft
+-!    !$           if((rhor(ifft,ispden)+0.0001_dp > zero).and.(lavnlr(ifft,ispden) >= zero)) then
+-!    !$              lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
+-!    !$           else
+-!    !$ !DEBUG
+-!    !$              !write(0,*) "DEBUG lavnl: for ifft=",ifft,"(rhor(ifft,ispden)+0.0001_dp)=",&
+-!    !$              !     & (rhor(ifft,ispden)+0.0001_dp),'lavnlr=',lavnlr(ifft,ispden)
+-!    !$ !ENDDEBUG
+-!    !$              lavnlr(ifft,ispden)=zero
+-!    !$              !lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
+-!    !$           end if
+-!    !$        end do
+-!    !$
+-!    !$        write(filapp,*) 'neww-lavnl-0.00'
+-!    !$        call out1line(filapp,lavnlr,dtset,0,0)
+-!    !$        write(filapp,*) 'neww-lavnl-0.25'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
+-!    !$        write(filapp,*) 'neww-lavnl-0.50'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
++!    !%        write(filapp,*) 'orig-lavnl-0.00'
++!    !%        call out1line(filapp,lavnlr,dtset,0,0)
++!    !%        write(filapp,*) 'orig-lavnl-0.25'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
++!    !%        write(filapp,*) 'orig-lavnl-0.50'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
++!    !%        lavnlr=lavnlr*(rhor+0.0001_dp)
++!    !% !DEBUG
++!    !%write(0,*) "DEBUG ================ lavnl ==========================="
++!    !%write(0,*) "DEBUG lavnl: density analysis"
++!    !% !ENDDEBUG
++!    !%        do ifft=1,dtset%nfft
++!    !%           if((rhor(ifft,ispden)+0.0001_dp > zero).and.(lavnlr(ifft,ispden) >= zero)) then
++!    !%              lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
++!    !%           else
++!    !% !DEBUG
++!    !%              !write(0,*) "DEBUG lavnl: for ifft=",ifft,"(rhor(ifft,ispden)+0.0001_dp)=",&
++!    !%              !     & (rhor(ifft,ispden)+0.0001_dp),'lavnlr=',lavnlr(ifft,ispden)
++!    !% !ENDDEBUG
++!    !%              lavnlr(ifft,ispden)=zero
++!    !%              !lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
++!    !%           end if
++!    !%        end do
++!    !%
++!    !%        write(filapp,*) 'neww-lavnl-0.00'
++!    !%        call out1line(filapp,lavnlr,dtset,0,0)
++!    !%        write(filapp,*) 'neww-lavnl-0.25'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
++!    !%        write(filapp,*) 'neww-lavnl-0.50'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
+ 
+ !    i1=1
+ !    i2=1
+diff -Naur abinit-6.4.2.bak/src/68_rsprc/lavnl.F90 abinit-6.4.2/src/68_rsprc/lavnl.F90
+--- src/68_rsprc/lavnl.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_rsprc/lavnl.F90	2011-01-05 09:30:40.000000000 +0000
+@@ -172,9 +172,9 @@
+ ! *************************************************************************
+ 
+ !***********************************************************************************
+-!Getting the localy averaged non-local potential                                ***
+-!$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
+-!**********************************************************************************
++! Getting the localy averaged non-local potential                                ***
++!%$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
++!***********************************************************************************
+ 
+  qphon=zero
+ 
+@@ -567,36 +567,36 @@
+      end do
+      lavnlr(:,ispden) = lavnlr(:,ispden)-min
+ !    lavnlr(:,ispden) = -min
+-!    !$        write(filapp,*) 'orig-lavnl-0.00'
+-!    !$        call out1line(filapp,lavnlr,dtset,0,0)
+-!    !$        write(filapp,*) 'orig-lavnl-0.25'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
+-!    !$        write(filapp,*) 'orig-lavnl-0.50'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
+-!    !$        lavnlr=lavnlr*(rhor+0.0001_dp)
+-!    !$ !DEBUG
+-!    !$write(0,*) "DEBUG ================ lavnl ==========================="
+-!    !$write(0,*) "DEBUG lavnl: density analysis"
+-!    !$ !ENDDEBUG
+-!    !$        do ifft=1,dtset%nfft
+-!    !$           if((rhor(ifft,ispden)+0.0001_dp > zero).and.(lavnlr(ifft,ispden) >= zero)) then
+-!    !$              lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
+-!    !$           else
+-!    !$ !DEBUG
+-!    !$              !write(0,*) "DEBUG lavnl: for ifft=",ifft,"(rhor(ifft,ispden)+0.0001_dp)=",&
+-!    !$              !     & (rhor(ifft,ispden)+0.0001_dp),'lavnlr=',lavnlr(ifft,ispden)
+-!    !$ !ENDDEBUG
+-!    !$              lavnlr(ifft,ispden)=zero
+-!    !$              !lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
+-!    !$           end if
+-!    !$        end do
+-!    !$
+-!    !$        write(filapp,*) 'neww-lavnl-0.00'
+-!    !$        call out1line(filapp,lavnlr,dtset,0,0)
+-!    !$        write(filapp,*) 'neww-lavnl-0.25'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
+-!    !$        write(filapp,*) 'neww-lavnl-0.50'
+-!    !$        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
++!    !%        write(filapp,*) 'orig-lavnl-0.00'
++!    !%        call out1line(filapp,lavnlr,dtset,0,0)
++!    !%        write(filapp,*) 'orig-lavnl-0.25'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
++!    !%        write(filapp,*) 'orig-lavnl-0.50'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
++!    !%        lavnlr=lavnlr*(rhor+0.0001_dp)
++!    !% !DEBUG
++!    !%write(0,*) "DEBUG ================ lavnl ==========================="
++!    !%write(0,*) "DEBUG lavnl: density analysis"
++!    !% !ENDDEBUG
++!    !%        do ifft=1,dtset%nfft
++!    !%           if((rhor(ifft,ispden)+0.0001_dp > zero).and.(lavnlr(ifft,ispden) >= zero)) then
++!    !%              lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
++!    !%           else
++!    !% !DEBUG
++!    !%              !write(0,*) "DEBUG lavnl: for ifft=",ifft,"(rhor(ifft,ispden)+0.0001_dp)=",&
++!    !%              !     & (rhor(ifft,ispden)+0.0001_dp),'lavnlr=',lavnlr(ifft,ispden)
++!    !% !ENDDEBUG
++!    !%              lavnlr(ifft,ispden)=zero
++!    !%              !lavnlr(ifft,ispden) = lavnlr(ifft,ispden)/(rhor(ifft,ispden)+0.0001_dp)
++!    !%           end if
++!    !%        end do
++!    !%
++!    !%        write(filapp,*) 'neww-lavnl-0.00'
++!    !%        call out1line(filapp,lavnlr,dtset,0,0)
++!    !%        write(filapp,*) 'neww-lavnl-0.25'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/4,dtset%ngfft(2)/4)
++!    !%        write(filapp,*) 'neww-lavnl-0.50'
++!    !%        call out1line(filapp,lavnlr,dtset,dtset%ngfft(1)/2,dtset%ngfft(2)/2)
+ 
+ !    i1=1
+ !    i2=1
+diff -Naur abinit-6.4.2.bak/src/68_rsprc/prctfvw1.F90 abinit-6.4.2/src/68_rsprc/prctfvw1.F90
+--- src/68_rsprc/prctfvw1.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_rsprc/prctfvw1.F90	2011-01-05 09:32:35.000000000 +0000
+@@ -209,9 +209,9 @@
+   type(cprj_type) :: cprj_dum(1,1)
+ ! *************************************************************************
+ !***********************************************************************************
+-!Getting the localy averaged non-local potential                                ***
+-!$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
+-!**********************************************************************************
++! Getting the localy averaged non-local potential                                ***
++!%$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
++!***********************************************************************************
+  znucl=1.0_dp
+  typat=1
+  xredcp = xred
+diff -Naur abinit-6.4.2.bak/src/68_rsprc/prctfvw2.F90 abinit-6.4.2/src/68_rsprc/prctfvw2.F90
+--- src/68_rsprc/prctfvw2.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_rsprc/prctfvw2.F90	2011-01-05 09:30:14.000000000 +0000
+@@ -228,8 +228,8 @@
+  xredcp = xred
+ !*************************************************************************
+ !***********************************************************************************
+-!Getting the localy averaged non-local potential                                ***
+-!$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
++! Getting the localy averaged non-local potential                                ***
++!%$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
+ !***********************************************************************************
+ !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
+  nfftotf=ngfftf(1)*ngfftf(2)*ngfftf(3)
+@@ -659,8 +659,8 @@
+ !******************************************************************
+ !Finding the density which minimizes the associated Energy      **
+ !******************************************************************
+-!!$  !compute the total charge number
+-!!$  !first switch to the sqrt of the density...
++!!%  !compute the total charge number
++!!%  !first switch to the sqrt of the density...
+  cplex=1;
+  option=1;
+  call dotprod_vn(cplex,& !complex density/pot
+@@ -676,46 +676,46 @@
+ &ucvol)          !cell volume
+ !enable the use of the functions eneofrho_tfw and deneofrho_tfw
+  Z=real(nint(Z),dp)
+-!!$  call init_eneofrho_tfw(dtset,dtset%intxc,dtset%ixc,psps%usepaw,n3xccc,ngfftf,&
+-!!$       &nfftf,nkxc,nspden,mpi_enreg,deltaW,gprimd,gsqcut,&
+-!!$       &lavnlfft,rhor,rprimd,ucvol,vout_unmixed,vpsp,vtrial,xccc3d,Z)
+-!!$  !minimizes Etfw with respect to sqrtrhor instead of rhor
+-!!$  !sqrtrhor(:,:)=two*Z/nfftf
+-!!$  !call random_number(rhor)
+-!!$  !eei=zero
+-!!$  !do ifft=1,nfft
+-!!$  !   eei=max(eei,sqrtrhor(ifft,1))
+-!!$  !end do
+-!!$  !eei=0.05_dp*eei
+-!!$  !call newdensity(eei,rhor,sqrtrhor)
+-!!$  call cgpr(eneofrho_tfw,deneofrho_tfw,newdensity,abs(deltae*real(0.001,dp)/etotal),55,sqrtrhor,dummy,dummy2)
+-!!$  !dummy=eneofrho_tfw(sqrtrhor)
+-!!$  ! free the dynamically allocated memory used by eneofrho_tfw and deneofrho_tfw
+-!!$  call end_eneofrho_tfw()
+-!!$  !new density from the minimised sqrtrhor
+-!!$  count=0
+-!!$  do ifft=1,nfftf
+-!!$     if (sqrtrhor(ifft,1)<zero) then
+-!!$        count=count+1
+-!!$     end if
+-!!$  end do
+-!!$  write(0,*) 'nombre sous zero',count
+-!!$  rhor(:,:)=sqrtrhor(:,:)*sqrtrhor(:,:)
+-!!$
+-!!$
+-!!$  write(0,*) 'n1 n2 n3',n1,n2,n3
+-!!$  i1=1
+-!!$  i2=1
+-!!$
+-!!$  do i3=1,n3
+-!!$  ifft=i1+n1*(i2-1+n2*(i3-1))
+-!!$  ifft2=i1+6+n1*(i2+6-1+n2*(i3-1))
+-!!$     write(79,*) deltaW(ifft,1),deltaW(ifft2,1),&
+-!!$     &         lavnlfft(ifft,1),lavnlfft(ifft2,1),&
+-!!$     &         vpsp(ifft)+vout_unmixed(ifft,1)-vtrial(ifft,1),&
+-!!$     &   vpsp(ifft2)+vout_unmixed(ifft2,1)-vtrial(ifft2,1),&
+-!!$     &         sqrtrhor(ifft,1),sqrtrhor(ifft2,1)
+-!!$  end do
++!!%  call init_eneofrho_tfw(dtset,dtset%intxc,dtset%ixc,psps%usepaw,n3xccc,ngfftf,&
++!!%       &nfftf,nkxc,nspden,mpi_enreg,deltaW,gprimd,gsqcut,&
++!!%       &lavnlfft,rhor,rprimd,ucvol,vout_unmixed,vpsp,vtrial,xccc3d,Z)
++!!%  !minimizes Etfw with respect to sqrtrhor instead of rhor
++!!%  !sqrtrhor(:,:)=two*Z/nfftf
++!!%  !call random_number(rhor)
++!!%  !eei=zero
++!!%  !do ifft=1,nfft
++!!%  !   eei=max(eei,sqrtrhor(ifft,1))
++!!%  !end do
++!!%  !eei=0.05_dp*eei
++!!%  !call newdensity(eei,rhor,sqrtrhor)
++!!%  call cgpr(eneofrho_tfw,deneofrho_tfw,newdensity,abs(deltae*real(0.001,dp)/etotal),55,sqrtrhor,dummy,dummy2)
++!!%  !dummy=eneofrho_tfw(sqrtrhor)
++!!%  ! free the dynamically allocated memory used by eneofrho_tfw and deneofrho_tfw
++!!%  call end_eneofrho_tfw()
++!!%  !new density from the minimised sqrtrhor
++!!%  count=0
++!!%  do ifft=1,nfftf
++!!%     if (sqrtrhor(ifft,1)<zero) then
++!!%        count=count+1
++!!%     end if
++!!%  end do
++!!%  write(0,*) 'nombre sous zero',count
++!!%  rhor(:,:)=sqrtrhor(:,:)*sqrtrhor(:,:)
++!!%
++!!%
++!!%  write(0,*) 'n1 n2 n3',n1,n2,n3
++!!%  i1=1
++!!%  i2=1
++!!%
++!!%  do i3=1,n3
++!!%  ifft=i1+n1*(i2-1+n2*(i3-1))
++!!%  ifft2=i1+6+n1*(i2+6-1+n2*(i3-1))
++!!%     write(79,*) deltaW(ifft,1),deltaW(ifft2,1),&
++!!%     &         lavnlfft(ifft,1),lavnlfft(ifft2,1),&
++!!%     &         vpsp(ifft)+vout_unmixed(ifft,1)-vtrial(ifft,1),&
++!!%     &   vpsp(ifft2)+vout_unmixed(ifft2,1)-vtrial(ifft2,1),&
++!!%     &         sqrtrhor(ifft,1),sqrtrhor(ifft2,1)
++!!%  end do
+ !******************************************************************
+ !Finding the density which minimizes the first term             **
+ !******************************************************************
+@@ -764,46 +764,46 @@
+ !vrespc=vrespc/diemix
+ !call laplacian(vrespc,lvres,ngfft,rprimd)
+ !lvres=-lvres*(two*two_pi)
+-!!$  write(0,*) 'efermi=',efermi
+-!!$  write(filapp,770) 'myoutput-opA',1
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,opA,xredcp,znucl)
+-!!$  write(filapp,770) 'myoutput-opB',1
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,opB,xredcp,znucl)
+-!!$  write(filapp,770) 'myoutput-partial',1
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,vpsp,xredcp,znucl)
+-!!$  write(filapp,770) 'myoutput-partial',2
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,vhartr,xredcp,znucl)
+-!!$  write(filapp,770) 'myoutput-partial',3
+-!!$  bla=seven/six*(alpha32*rhor(:,:))**(two_thirds)
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,bla,xredcp,znucl)
+-!!$  write(filapp,770) 'myoutput-partial',4
+-!!$  bla(:,1)=2*kxc(:,1)*rhor(:,1)
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,bla,xredcp,znucl)
+-!!$  write(filapp,770) 'myoutput-partial',5
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,lavnlfft,xredcp,znucl)
+-!!$  write(filapp,770) 'myoutput-partial',6
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,vxc,xredcp,znucl)
+-!!$bla=efermi
+-!!$  write(filapp,770) 'myoutput-partial',7
+-!!$  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
+-!!$       &  rhor,rprimd,&
+-!!$       &  typat,ucvol,bla,xredcp,znucl)
++!!%  write(0,*) 'efermi=',efermi
++!!%  write(filapp,770) 'myoutput-opA',1
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,opA,xredcp,znucl)
++!!%  write(filapp,770) 'myoutput-opB',1
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,opB,xredcp,znucl)
++!!%  write(filapp,770) 'myoutput-partial',1
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,vpsp,xredcp,znucl)
++!!%  write(filapp,770) 'myoutput-partial',2
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,vhartr,xredcp,znucl)
++!!%  write(filapp,770) 'myoutput-partial',3
++!!%  bla=seven/six*(alpha32*rhor(:,:))**(two_thirds)
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,bla,xredcp,znucl)
++!!%  write(filapp,770) 'myoutput-partial',4
++!!%  bla(:,1)=2*kxc(:,1)*rhor(:,1)
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,bla,xredcp,znucl)
++!!%  write(filapp,770) 'myoutput-partial',5
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,lavnlfft,xredcp,znucl)
++!!%  write(filapp,770) 'myoutput-partial',6
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,vxc,xredcp,znucl)
++!!%bla=efermi
++!!%  write(filapp,770) 'myoutput-partial',7
++!!%  call out1dm(filapp,natom,nfftf,ngfftf,nspden,ntypat,&
++!!%       &  rhor,rprimd,&
++!!%       &  typat,ucvol,bla,xredcp,znucl)
+  a=1.e9
+  b=1.e9
+  do ifft=1,nfftf
+@@ -837,17 +837,17 @@
+ !******************************************************************
+  cplex=1;
+  option=1;
+-!!$  call dotprod_vn(cplex,& !complex density/pot
+-!!$       &phi2,&           !density
+-!!$       &a,&              !resulting dorproduct integrated over r
+-!!$       &doti,&           !imaginary part of the integral
+-!!$       &mpi_enreg,&      !
+-!!$       &size(rhor,1),&   !number of localy(cpu) attributed grid point
+-!!$       &nfftotf,&        !real total number of grid point
+-!!$       &size(rhor,2),&   !nspden
+-!!$       &option,&         !1=compute only the real part 2=compute also the imaginary part
+-!!$       &phi2,&           !density
+-!!$       &ucvol)           !cell volume
++!!%  call dotprod_vn(cplex,& !complex density/pot
++!!%       &phi2,&           !density
++!!%       &a,&              !resulting dorproduct integrated over r
++!!%       &doti,&           !imaginary part of the integral
++!!%       &mpi_enreg,&      !
++!!%       &size(rhor,1),&   !number of localy(cpu) attributed grid point
++!!%       &nfftotf,&        !real total number of grid point
++!!%       &size(rhor,2),&   !nspden
++!!%       &option,&         !1=compute only the real part 2=compute also the imaginary part
++!!%       &phi2,&           !density
++!!%       &ucvol)           !cell volume
+  call dotprod_vn(cplex,& !complex density/pot
+ &phi2,&           !density
+ &b,&              !resulting dorproduct integrated over r
+@@ -870,19 +870,19 @@
+ &option,&         !1=compute only the real part 2=compute also the imaginary part
+ &sqrtrhor,&           !density
+ &ucvol)           !cell volume
+-!!$  discr= b*b-four*a*c
+-!!$  if(discr > zero) then
+-!!$     discr=sqrt(discr)
+-!!$     fermi1p=(-b+discr)/(two*a)
+-!!$     fermi1m=(-b-discr)/(two*a)
+-!!$  else if (discr == zero) then
+-!!$     fermi1p=-b/(two*a)
+-!!$     fermi1m=-b/(two*a)
+-!!$  else
+-!!$     fermi1p=zero
+-!!$     fermi1m=zero
+-!!$  end if
+-!!$  dfermi=merge(fermi1p,fermi1m,abs(fermi1p) <= abs(fermi1m))
++!!%  discr= b*b-four*a*c
++!!%  if(discr > zero) then
++!!%     discr=sqrt(discr)
++!!%     fermi1p=(-b+discr)/(two*a)
++!!%     fermi1m=(-b-discr)/(two*a)
++!!%  else if (discr == zero) then
++!!%     fermi1p=-b/(two*a)
++!!%     fermi1m=-b/(two*a)
++!!%  else
++!!%     fermi1p=zero
++!!%     fermi1m=zero
++!!%  end if
++!!%  dfermi=merge(fermi1p,fermi1m,abs(fermi1p) <= abs(fermi1m))
+  dfermi = -c/b !since the term in a must be neglected !(??)
+  sqrtrhor(:,:)=phi1+dfermi*phi2
+ !****************************************************************** !
+@@ -1056,115 +1056,115 @@
+ !vtrial fourier
+ !vnewfourier
+ !ref
+-!!$  count=0
+-!!$  count1=0
+-!!$  count2=0
+-!!$  count3=0
+-!!$  count4=0
+-!!$  count5=0
+-!!$  count6=0
+-!!$  count7=0
+-!!$  count8=0
+-!!$  count9=0
+-!!$  count10=0
+-!!$  count11=0
+-!!$  count12=0
+-!!$  count13=0
+-!!$  count14=0
+-!!$  count15=0
+-!!$  count16=0
+-!!$  count17=0
+-!!$  count18=0
+-!!$  do ispden=1,nspden
+-!!$     do ifft=1,nfftf
+-!!$        aa=vrefg(:,ifft,ispden)
+-!!$        bb=vin_oldfourier(:,ifft,ispden)-aa
+-!!$        cc=newvoutfourier(:,ifft,ispden)-aa
+-!!$        principal: if((g2cart(ifft).lt.one).and.(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-5) then
+-!!$           ong: if(g2cart(ifft).lt.one*0.2) then
+-!!$              onv2:if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.2d-1) then
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count15=count15+1
+-!!$                 else
+-!!$                    count16=count16+1
+-!!$                 end if
+-!!$              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-2) then
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count1=count1+1
+-!!$                 else
+-!!$                    count2=count2+1
+-!!$                 end if
+-!!$
+-!!$              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-3) then
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count3=count3+1
+-!!$                 else
+-!!$                    count4=count4+1
+-!!$                 end if
+-!!$              else
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count5=count5+1
+-!!$                 else
+-!!$                    count6=count6+1
+-!!$                 end if
+-!!$              end if onv2
+-!!$           else
+-!!$              onv22:if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.2d-1) then
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count17=count17+1
+-!!$                 else
+-!!$                    count18=count18+1
+-!!$                 end if
+-!!$              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-2) then
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count1=count1+1
+-!!$                 else
+-!!$                    count2=count2+1
+-!!$                 end if
+-!!$              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-3) then
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count9=count9+1
+-!!$                 else
+-!!$                    count10=count10+1
+-!!$                 end if
+-!!$              else
+-!!$                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
+-!!$                    count11=count11+1
+-!!$                 else
+-!!$                    count12=count12+1
+-!!$                 end if
+-!!$              end if onv22
+-!!$           end if ong
+-!!$        end if principal
+-!!$        if(abs(vref(ifft,ispden)-vtrial(ifft,ispden)).gt.&
+-!!$             &abs(vref(ifft,ispden)-vtrialold(ifft,ispden))) then
+-!!$           count14=count14+1
+-!!$        else
+-!!$           count13=count13+1
+-!!$        end if
+-!!$
+-!!$     end do
+-!!$  end do
+-!!$  write(1111,*),'---------------------------------------------------------------------------------'
+-!!$  write(1111,*),'|           g2<0.2          |           g2<1.0          |'
+-!!$  write(1111,*),'|     ok      |    bad      |     ok      |    bad      |  '
+-!!$  write(1111,*),'|',count16,'|',count15,'|',count18,'|',count17,'|  ','||V||>2e-1'
+-!!$  write(1111,*),'|',count2,'|',count1,'|',count8,'|',count7,'|  ','||V||>5e-2'
+-!!$  write(1111,*),'|',count4,'|',count3,'|',count10,'|',count9,'|  ','||V||>5e-3'
+-!!$  write(1111,*),'|',count6,'|',count5,'|',count12,'|',count11,'|  ','||V||>5e-5'
+-!!$  write(1111,*),'---------------------------------------------------------------------------------'
+-!!$  write(0,*),'---------------------------------------------------------------------------------'
+-!!$  write(0,*),'|           g2<0.2          |           g2<1.0          |'
+-!!$  write(0,*),'|     ok      |    bad      |     ok      |    bad      |  '
+-!!$  write(0,*),'|',count16,'|',count15,'|',count18,'|',count17,'|  ','||V||>2e-1'
+-!!$  write(0,*),'|',count2,'|',count1,'|',count8,'|',count7,'|  ','||V||>5e-2'
+-!!$  write(0,*),'|',count4,'|',count3,'|',count10,'|',count9,'|  ','||V||>5e-3'
+-!!$  write(0,*),'|',count6,'|',count5,'|',count12,'|',count11,'|  ','||V||>5e-5'
+-!!$  write(0,*),'---------------------------------------------------------------------------------'
+-!!$
+-!!$
+-!!$  write(0,*) 'les mauvais r:',count14,'les bons r',count13
+-!!$  write(1111,*) 'les mauvais r:',count14,'les bons r',count13
++!!%  count=0
++!!%  count1=0
++!!%  count2=0
++!!%  count3=0
++!!%  count4=0
++!!%  count5=0
++!!%  count6=0
++!!%  count7=0
++!!%  count8=0
++!!%  count9=0
++!!%  count10=0
++!!%  count11=0
++!!%  count12=0
++!!%  count13=0
++!!%  count14=0
++!!%  count15=0
++!!%  count16=0
++!!%  count17=0
++!!%  count18=0
++!!%  do ispden=1,nspden
++!!%     do ifft=1,nfftf
++!!%        aa=vrefg(:,ifft,ispden)
++!!%        bb=vin_oldfourier(:,ifft,ispden)-aa
++!!%        cc=newvoutfourier(:,ifft,ispden)-aa
++!!%        principal: if((g2cart(ifft).lt.one).and.(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-5) then
++!!%           ong: if(g2cart(ifft).lt.one*0.2) then
++!!%              onv2:if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.2d-1) then
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count15=count15+1
++!!%                 else
++!!%                    count16=count16+1
++!!%                 end if
++!!%              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-2) then
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count1=count1+1
++!!%                 else
++!!%                    count2=count2+1
++!!%                 end if
++!!%
++!!%              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-3) then
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count3=count3+1
++!!%                 else
++!!%                    count4=count4+1
++!!%                 end if
++!!%              else
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count5=count5+1
++!!%                 else
++!!%                    count6=count6+1
++!!%                 end if
++!!%              end if onv2
++!!%           else
++!!%              onv22:if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.2d-1) then
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count17=count17+1
++!!%                 else
++!!%                    count18=count18+1
++!!%                 end if
++!!%              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-2) then
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count1=count1+1
++!!%                 else
++!!%                    count2=count2+1
++!!%                 end if
++!!%              else if(sqrt(vrefg(1,ifft,ispden)**2+vrefg(2,ifft,ispden)).gt.5d-3) then
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count9=count9+1
++!!%                 else
++!!%                    count10=count10+1
++!!%                 end if
++!!%              else
++!!%                 if(bb(1)**2+bb(2)**2.lt.cc(1)**2+cc(2)**2) then
++!!%                    count11=count11+1
++!!%                 else
++!!%                    count12=count12+1
++!!%                 end if
++!!%              end if onv22
++!!%           end if ong
++!!%        end if principal
++!!%        if(abs(vref(ifft,ispden)-vtrial(ifft,ispden)).gt.&
++!!%             &abs(vref(ifft,ispden)-vtrialold(ifft,ispden))) then
++!!%           count14=count14+1
++!!%        else
++!!%           count13=count13+1
++!!%        end if
++!!%
++!!%     end do
++!!%  end do
++!!%  write(1111,*),'---------------------------------------------------------------------------------'
++!!%  write(1111,*),'|           g2<0.2          |           g2<1.0          |'
++!!%  write(1111,*),'|     ok      |    bad      |     ok      |    bad      |  '
++!!%  write(1111,*),'|',count16,'|',count15,'|',count18,'|',count17,'|  ','||V||>2e-1'
++!!%  write(1111,*),'|',count2,'|',count1,'|',count8,'|',count7,'|  ','||V||>5e-2'
++!!%  write(1111,*),'|',count4,'|',count3,'|',count10,'|',count9,'|  ','||V||>5e-3'
++!!%  write(1111,*),'|',count6,'|',count5,'|',count12,'|',count11,'|  ','||V||>5e-5'
++!!%  write(1111,*),'---------------------------------------------------------------------------------'
++!!%  write(0,*),'---------------------------------------------------------------------------------'
++!!%  write(0,*),'|           g2<0.2          |           g2<1.0          |'
++!!%  write(0,*),'|     ok      |    bad      |     ok      |    bad      |  '
++!!%  write(0,*),'|',count16,'|',count15,'|',count18,'|',count17,'|  ','||V||>2e-1'
++!!%  write(0,*),'|',count2,'|',count1,'|',count8,'|',count7,'|  ','||V||>5e-2'
++!!%  write(0,*),'|',count4,'|',count3,'|',count10,'|',count9,'|  ','||V||>5e-3'
++!!%  write(0,*),'|',count6,'|',count5,'|',count12,'|',count11,'|  ','||V||>5e-5'
++!!%  write(0,*),'---------------------------------------------------------------------------------'
++!!%
++!!%
++!!%  write(0,*) 'les mauvais r:',count14,'les bons r',count13
++!!%  write(1111,*) 'les mauvais r:',count14,'les bons r',count13
+ !stop
+ !----------------------------------------------------------------------------------------------------------
+ !last FREE
+diff -Naur abinit-6.4.2.bak/src/68_rsprc/prctfw3.F90 abinit-6.4.2/src/68_rsprc/prctfw3.F90
+--- src/68_rsprc/prctfw3.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/68_rsprc/prctfw3.F90	2011-01-05 09:33:04.000000000 +0000
+@@ -170,9 +170,9 @@
+ ! *************************************************************************
+ 
+ !***********************************************************************************
+-!Getting the localy averaged non-local potential                                ***
+-!$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
+-!**********************************************************************************
++! Getting the localy averaged non-local potential                                ***
++!%$Vnl(r) = [\sum_{n,k} f_{n,k} \psi_{n,k}(r) (Vnl(r,r') |\psi_{n,k}(r')>)]/n(r)$***
++!***********************************************************************************
+  rhor=rhor_in
+  qphon=zero
+  xred2=xred
+diff -Naur abinit-6.4.2.bak/src/69_bse/cexch_haydock.F90 abinit-6.4.2/src/69_bse/cexch_haydock.F90
+--- src/69_bse/cexch_haydock.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/69_bse/cexch_haydock.F90	2011-01-05 16:47:03.000000000 +0000
+@@ -185,7 +185,7 @@
+ ! There is no need to shift the G-sphere as we only have vertical transitions.
+  mgfft_osc   = MAXVAL(ngfft_osc(1:3))
+  fftalga_osc = ngfft_osc(7)/100 !; fftalgc_osc=MOD(ngfft_osc(7),10)
+- use_padfft=0; !$if (fftalga_osc==3) use_padfft=1  ! Padded FFTW3 is safe instead!
++ use_padfft=0; !%if (fftalga_osc==3) use_padfft=1  ! Padded FFTW3 is safe instead!
+  allocate(gbound(2*mgfft_osc+8,2*use_padfft))
+  if (use_padfft==1) call sphereboundary(gbound,1,Gsph_Max%gvec,mgfft_osc,BSp%npweps)
+ 
+@@ -315,8 +315,8 @@
+    do ic=BSp%lomo,BSp%nbnds ! loop over band C 
+ !   Symmetry index for IBZ = 1 (in ktabr(:,1)
+ 
+-    !$call wfd_get_ur(Wfd,iv,ik,spin,wfr1)
+-    !$call wfd_get_ur(Wfd,ic,ik,spin,wfr2)
++    !%call wfd_get_ur(Wfd,iv,ik,spin,wfr1)
++    !%call wfd_get_ur(Wfd,ic,ik,spin,wfr2)
+ 
+     call rho_tw_g(paral_kgb,nspinor,BSp%npweps,nfftot_osc,ngfft_osc,map2sphere,use_padfft,igfftg0,gbound,&
+ &     wfnr(:,iv,ik),1,ktabr(:,1),ph_mkt,spinrot_k,& 
+@@ -533,8 +533,8 @@
+       do iv=BSp%lomo,BSp%homo    ! loop over band V
+         do ivp=BSp%lomo,BSp%homo  ! loop over band VP
+ 
+-          !$call wfd_get_ur(Wfd,ivp,ikp_ibz,spin,wfr1)
+-          !$call wfd_get_ur(Wfd,iv ,ik_ibz ,spin,wfr2)
++          !%call wfd_get_ur(Wfd,ivp,ikp_ibz,spin,wfr1)
++          !%call wfd_get_ur(Wfd,iv ,ik_ibz ,spin,wfr2)
+ 
+           call rho_tw_g(paral_kgb,nspinor,BSp%npweps,nfftot_osc,ngfft_osc,map2sphere,use_padfft,igfftg0,gbound,&
+ &           wfnr(:,ivp,ikp_ibz),itim_kp,ktabr_kp,ph_mkpt,spinrot_kp,&
+@@ -556,8 +556,8 @@
+       do icp=BSp%lumo,BSp%nbnds  ! loop over band C_prime
+         do ic=BSp%lumo,BSp%nbnds  ! loop over band C
+ 
+-          !$call wfd_get_ur(Wfd,icp,ikp_ibz,spin,wfr1)
+-          !$call wfd_get_ur(Wfd,ic ,ik_ibz ,spin,wfr2)
++          !%call wfd_get_ur(Wfd,icp,ikp_ibz,spin,wfr1)
++          !%call wfd_get_ur(Wfd,ic ,ik_ibz ,spin,wfr2)
+ 
+           call rho_tw_g(paral_kgb,nspinor,BSp%npweps,nfftot_osc,ngfft_osc,map2sphere,use_padfft,igfftg0,gbound,&
+ &           wfnr(:,icp,ikp_ibz),itim_kp,ktabr_kp,ph_mkpt,spinrot_kp,&
+@@ -613,7 +613,7 @@
+ 
+ !           then sum on G': sum_G' (rho_c'c(G) W*(G,G'))* rho_v'v(G')
+             http = -faq * DOT_PRODUCT(ctccp(:),rhxtwg_vv(:,ivp,iv))
+-            !$http = -faq * xdotc(BSp%npweps,ctccp,1,rhxtwg_vv(:,ivp,iv),1)
++            !%http = -faq * xdotc(BSp%npweps,ctccp,1,rhxtwg_vv(:,ivp,iv),1)
+ 
+             bsh_k(irv,irc,ikp,irvp,ircp) = http
+ !           write(*,'("ik,ikp,iv,ivp,ic,icp, it, itp",8I3)')ik,ikp,iv,ivp,ic,icp, it, itp
+@@ -663,7 +663,7 @@
+ 
+        ! sum over G
+        ctemp = DOT_PRODUCT(rhotwg1(2:),rhotwg2(2:))
+-       !$ctemp = XDOTC(BSp%npweps-1,rhotwg1(2:),1,rhotwg2(2:),1)
++       !%ctemp = XDOTC(BSp%npweps-1,rhotwg1(2:),1,rhotwg2(2:),1)
+        ctemp = two * ctemp * faq
+                
+        bsh_k(irv,irc,ikp,irvp,ircp)=bsh_k(irv,irc,ikp,irvp,ircp)+ctemp
+diff -Naur abinit-6.4.2.bak/src/69_bse/exceig.F90 abinit-6.4.2/src/69_bse/exceig.F90
+--- src/69_bse/exceig.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/69_bse/exceig.F90	2011-01-05 16:48:51.000000000 +0000
+@@ -155,11 +155,11 @@
+ 
+    ! Check the size of hamiltonian matrix
+    ! Not possible anymore, this info should be reported in the header.
+-   !$read(hreso_unt,rec=((nh*nh + nh)/2)+1) itemp
+-   !$if(itemp/=nh) then
+-   !$ write(msg,'(a,i0,a,i0,a)')'nt should be ',nh,' but read ',itemp,' wrong file fort.55 *.exh'
+-   !$ MSG_ERROR(msg)
+-   !$end if
++   !%read(hreso_unt,rec=((nh*nh + nh)/2)+1) itemp
++   !%if(itemp/=nh) then
++   !% write(msg,'(a,i0,a,i0,a)')'nt should be ',nh,' but read ',itemp,' wrong file fort.55 *.exh'
++   !% MSG_ERROR(msg)
++   !%end if
+ 
+    ! Construct full excitonic Hamiltonian using Hermiticity. file is always in double precision.
+    do itp=1,nh
+@@ -181,7 +181,7 @@
+      call wrtout(std_out," Partial diagonalization via XHEEVX","COLL")
+      abstol=zero; ldz=nh
+      allocate(exc_vec(ldz,nh)) ! TODO Single precision is not available.
+-     !$call xheevx("Vectors","All","Upper",nh,exc_mat,vl,vu,il,iu,abstol,mene_found,exc_ene,exc_vec,ldz)
++     !%call xheevx("Vectors","All","Upper",nh,exc_mat,vl,vu,il,iu,abstol,mene_found,exc_ene,exc_vec,ldz)
+      exc_mat = exc_vec
+      deallocate(exc_vec)
+    end if
+@@ -287,28 +287,28 @@
+ 
+    call wrtout(std_out,' Writing eigenvalues/vectors on file: '//TRIM(filbseig),"COLL")
+ 
+-   !$ Very inefficient coding. Allocate global array to store the distributed eigenvectors  
+-   !$ Fill the matrix, then the master node writes the final results on file 
+-   !$allocate(exc_mat_dpc(nh,nh), STAT=istat)
+-   !$ABI_CHECK(istat==0,'out of memory: excitonic eigenvectors')
+-   !$exc_mat_dpc = czero
+-
+-   !$call slk_matrix_to_global_dpc_2D(Slk_vec,"All",exc_mat_dpc)
+-   !$call xsum_master(exc_mat_dpc,master,spaceComm,ierr)
+-
+-   !$if (my_rank==master) then
+-   !$ eig_unt=get_unit()
+-   !$ open(eig_unt,file=filbseig,form='unformatted')
+-   !$                                                      
+-   !$ write(eig_unt) nh
+-   !$ write(eig_unt) CMPLX(exc_ene(1:nh),kind=dpc)
+-   !$ do mi=1,nh
+-   !$  write(eig_unt) CMPLX(exc_mat_dpc(1:nh,mi),kind=dpc)
+-   !$ end do
+-   !$ close(eig_unt)
+-   !$end if
+-   !$
+-   !$deallocate(exc_mat_dpc)
++   !% Very inefficient coding. Allocate global array to store the distributed eigenvectors  
++   !% Fill the matrix, then the master node writes the final results on file 
++   !%allocate(exc_mat_dpc(nh,nh), STAT=istat)
++   !%ABI_CHECK(istat==0,'out of memory: excitonic eigenvectors')
++   !%exc_mat_dpc = czero
++
++   !%call slk_matrix_to_global_dpc_2D(Slk_vec,"All",exc_mat_dpc)
++   !%call xsum_master(exc_mat_dpc,master,spaceComm,ierr)
++
++   !%if (my_rank==master) then
++   !% eig_unt=get_unit()
++   !% open(eig_unt,file=filbseig,form='unformatted')
++   !%                                                      
++   !% write(eig_unt) nh
++   !% write(eig_unt) CMPLX(exc_ene(1:nh),kind=dpc)
++   !% do mi=1,nh
++   !%  write(eig_unt) CMPLX(exc_mat_dpc(1:nh,mi),kind=dpc)
++   !% end do
++   !% close(eig_unt)
++   !%end if
++   !%
++   !%deallocate(exc_mat_dpc)
+ 
+    ! Write distributed matrix on file tmp_fname using streams instead of Fortran records.
+    tmp_fname=pick_aname()
+@@ -661,8 +661,8 @@
+ 
+  ! Check the size of hamiltonian matrix
+  ! Not possible anymore, this info should be reported in the header.
+- !$read(hreso_unt,rec=((nh*nh + nh)/2)+1) itemp
+- !$ABI_CHECK(itemp==nh,'Wrong resonant file')
++ !%read(hreso_unt,rec=((nh*nh + nh)/2)+1) itemp
++ !%ABI_CHECK(itemp==nh,'Wrong resonant file')
+ 
+  do itp=1,nh
+    do it=1,itp
+@@ -687,8 +687,8 @@
+ 
+  ! Check the size of hamiltonian matrix
+  ! Not possible anymore, this info should be reported in the header.
+- !$read(hcoup_unt,rec=((nh*nh + nh)/2)+1) itemp
+- !$ABI_CHECK(itemp==nh,'wrong file .exc')
++ !%read(hcoup_unt,rec=((nh*nh + nh)/2)+1) itemp
++ !%ABI_CHECK(itemp==nh,'wrong file .exc')
+ 
+  do itp=1,nh
+    do it=1,itp
+@@ -970,28 +970,28 @@
+   ! Fill the matrix, then master node writes the final results on file 
+   ! (better coding requires either scaLAPACK tool or MPI-IO (the later should be much faster).
+ 
+-  !$allocate(ovlp_dpc(nh2,nh2), STAT=istat)
+-  !$ABI_CHECK(istat==0,'out of memory ovlp_dpc')
+-  !$ovlp_dpc = czero
+-
+-  !$call slk_matrix_to_global_dpc_2D(Slk_mat,"All",ovlp_dpc)
+-  !$call xsum_master(ovlp_dpc,master,spaceComm,ierr)
+-
+-  !$if (my_rank==master) then
+-  !$ msg=' Writing overlap matrix s^-1 on file '//TRIM(BS_files%exovl)
+-  !$ call wrtout(std_out,msg,"PERS")
+-
+-  !$ ovlp_unt = get_unit()
+-  !$ open(unit=ovlp_unt,file=BS_files%exovl,form='unformatted',iostat=ios)
+-
+-  !$ write(ovlp_unt) nh2
+-  !$ do mi=1,nh2 ! write overlap matrix s^-1(mi,:)
+-  !$  write(ovlp_unt) CMPLX(ovlp_dpc(mi,:),kind=dpc)
+-  !$ end do
+-  !$ close(ovlp_unt)
+-  !$end if
++  !%allocate(ovlp_dpc(nh2,nh2), STAT=istat)
++  !%ABI_CHECK(istat==0,'out of memory ovlp_dpc')
++  !%ovlp_dpc = czero
++
++  !%call slk_matrix_to_global_dpc_2D(Slk_mat,"All",ovlp_dpc)
++  !%call xsum_master(ovlp_dpc,master,spaceComm,ierr)
++
++  !%if (my_rank==master) then
++  !% msg=' Writing overlap matrix s^-1 on file '//TRIM(BS_files%exovl)
++  !% call wrtout(std_out,msg,"PERS")
++
++  !% ovlp_unt = get_unit()
++  !% open(unit=ovlp_unt,file=BS_files%exovl,form='unformatted',iostat=ios)
++
++  !% write(ovlp_unt) nh2
++  !% do mi=1,nh2 ! write overlap matrix s^-1(mi,:)
++  !%  write(ovlp_unt) CMPLX(ovlp_dpc(mi,:),kind=dpc)
++  !% end do
++  !% close(ovlp_unt)
++  !%end if
+ 
+-  !$deallocate(ovlp_dpc)
++  !%deallocate(ovlp_dpc)
+ 
+   ! Use MPI-IO
+   ! Write distributed matrix on file tmp_fname using streams instead of Fortran records.
+diff -Naur abinit-6.4.2.bak/src/69_bse/setup_bse.F90 abinit-6.4.2/src/69_bse/setup_bse.F90
+--- src/69_bse/setup_bse.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/69_bse/setup_bse.F90	2011-01-05 16:50:53.000000000 +0000
+@@ -362,7 +362,7 @@
+  CASE (2)
+    msg= " Model dielectric function not yet coded"
+    MSG_ERROR(msg)
+-   !$BSp%wtype = 'wmodel' 
++   !%BSp%wtype = 'wmodel' 
+  CASE DEFAULT
+    write(msg,'(a,i0)')" Wrong bs_coulomb_term: ",Dtset%bs_coulomb_term
+    MSG_ERROR(msg)
+@@ -715,7 +715,7 @@
+ 
+  !TODO call update_occ here
+  ! Occupancies might be zero if NSCF 
+- !$call update_occ(KS_BSt,fixmom,stmbias,Dtset%prtvol)
++ !%call update_occ(KS_BSt,fixmom,stmbias,Dtset%prtvol)
+ 
+  call print_bandstructure(KS_BSt,"Band structure read from the KSS file",unit=std_out,prtvol=Dtset%prtvol)
+ 
+diff -Naur abinit-6.4.2.bak/src/77_ddb/eliashberg_1d.F90 abinit-6.4.2/src/77_ddb/eliashberg_1d.F90
+--- src/77_ddb/eliashberg_1d.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/77_ddb/eliashberg_1d.F90	2011-01-05 09:36:08.000000000 +0000
+@@ -108,10 +108,10 @@
+ 
+ !
+ !1) use linearized Eliashberg equation to find Tc
+-!$ \sum_j \mathbf{M}_{ij} \Delta_j = \zeta \cdot \Delta_i $ $i,j = 1 .. n_{\mathrm{Matsubara}}$
+-!$\zeta = 1$ gives T$_c$ $\beta = \frac{1}{\mathrm{T}}$ $\omega_i = (2 i + 1) \pi \mathrm{T}$
+-!$ \mathbf{M}_{ij} = \frac{\pi}{\beta} \frac{\lambda (\omega_i - \omega_j)}{Z (\omega_i)}$
+-!$ Z (\omega_i) = 1 + \frac{\pi}{\beta \omega_i} \sum_j \lambda(\omega_i - \omega_j) \mathrm{sgn}(\omega_j)$
++!%$ \sum_j \mathbf{M}_{ij} \Delta_j = \zeta \cdot \Delta_i $ $i,j = 1 .. n_{\mathrm{Matsubara}}$
++!%$\zeta = 1$ gives T$_c$ $\beta = \frac{1}{\mathrm{T}}$ $\omega_i = (2 i + 1) \pi \mathrm{T}$
++!%$ \mathbf{M}_{ij} = \frac{\pi}{\beta} \frac{\lambda (\omega_i - \omega_j)}{Z (\omega_i)}$
++!%$ Z (\omega_i) = 1 + \frac{\pi}{\beta \omega_i} \sum_j \lambda(\omega_i - \omega_j) \mathrm{sgn}(\omega_j)$
+ !
+ 
+ !initial guess for T$_c$ in Hartree (1Ha =3.067e5 K)
+diff -Naur abinit-6.4.2.bak/src/77_ddb/m_eph.F90 abinit-6.4.2/src/77_ddb/m_eph.F90
+--- src/77_ddb/m_eph.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/77_ddb/m_eph.F90	2011-01-05 09:34:59.000000000 +0000
+@@ -233,10 +233,10 @@
+  public :: nullify_fermi_surface
+  public :: destroy_fermi_surface
+ 
+- !$public :: init_fermi_surface
+- !$public :: wannier_interpolate_fermi_surface
+- !$public :: get_fs_ibz
+- !$public :: bxsf_write_fermi_surface
++ !%public :: init_fermi_surface
++ !%public :: wannier_interpolate_fermi_surface
++ !%public :: get_fs_ibz
++ !%public :: bxsf_write_fermi_surface
+ 
+  ! example:
+  !type(fermi_surface_type),allocatable :: Fsurf(:)
+@@ -325,9 +325,9 @@
+ ! Bound Methods: 
+  public :: nullify_gkk
+  public :: destroy_gkk
+- !$ init_gkk
+- !$ read_gkk_from_file
+- !$ get_gkk_full_fsbz        ! complete gkk on the full FS BZ.
++ !% init_gkk
++ !% read_gkk_from_file
++ !% get_gkk_full_fsbz        ! complete gkk on the full FS BZ.
+ 
+  interface nullify_gkk
+   module procedure nullify_gkk_0D
+@@ -379,9 +379,9 @@
+ ! Bound Methods: 
+  public :: nullify_gkk_handler
+  public :: destroy_gkk_handler
+- !$init_gkk_handler(Gkk,FSurf,Cryst,Cryst,qpt,fname)
+- !$get_gammaq
+- !$symmetrize_gkk_over_perts
++ !%init_gkk_handler(Gkk,FSurf,Cryst,Cryst,qpt,fname)
++ !%get_gammaq
++ !%symmetrize_gkk_over_perts
+ 
+  interface nullify_gkk_handler
+   module procedure nullify_gkk_handler_0D
+diff -Naur abinit-6.4.2.bak/src/77_suscep/prtsusd.F90 abinit-6.4.2/src/77_suscep/prtsusd.F90
+--- src/77_suscep/prtsusd.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/77_suscep/prtsusd.F90	2011-01-05 08:58:59.000000000 +0000
+@@ -86,7 +86,7 @@
+  allocate(sus_poly(npw_tiny))
+ 
+ !Scaling for e-e interaction energy
+-!$-1/2\pi$ [fluctuation-dissipation theorem] * $4\pi$ [Coulomb interaction] = -2
++!% $-1/2\pi$ [fluctuation-dissipation theorem] * $4\pi$ [Coulomb interaction] = -2
+  scalefactor=-2._dp ! -1/2pi[fluctuation-dissipation theorem] * 4pi[Coulomb interaction]
+ 
+ !Directional extrapolation followed by spatial average for G=0 term, using polynomial
+diff -Naur abinit-6.4.2.bak/src/95_drive/bethe_salpeter.F90 abinit-6.4.2/src/95_drive/bethe_salpeter.F90
+--- src/95_drive/bethe_salpeter.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/95_drive/bethe_salpeter.F90	2011-01-05 09:03:45.000000000 +0000
+@@ -567,7 +567,7 @@
+ 
+ !TODO this has to be done in a better way, moreover wont work for PAW
+ !Check Vcp!
+-!$ call cutoff_density(ngfftf,Dtset%nspden,Dtset%nsppol,Vcp,ks_rhor,MPI_enreg)
++!% call cutoff_density(ngfftf,Dtset%nspden,Dtset%nsppol,Vcp,ks_rhor,MPI_enreg)
+ !
+ !=== Additional computation for PAW ===
+  if (Dtset%usepaw==1) then
+diff -Naur abinit-6.4.2.bak/src/95_drive/gw_driver.F90 abinit-6.4.2/src/95_drive/gw_driver.F90
+--- src/95_drive/gw_driver.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/95_drive/gw_driver.F90	2011-01-05 09:22:30.000000000 +0000
+@@ -171,7 +171,7 @@
+      if (read_qps1/=0) then ! Change the default.
+        Dtfil%fnameabi_qps=filqps_it1
+      else                   ! Save the IT1 name for GW0 runs
+-       !$filqps_it1=TRIM(filnam_ds(3))//'_QPS'
++       !%filqps_it1=TRIM(filnam_ds(3))//'_QPS'
+        filqps_it1=Dtfil%fnameabi_qps
+      end if
+      !    
+@@ -182,7 +182,7 @@
+      if (read_scr1/=0) then  ! Change the default.
+        Dtfil%fnameabi_scr=filscr_it1
+      else                    ! Save the IT1 name for GW0 runs
+-       !$filscr_it1=TRIM(filnam_ds(3))//'_SCR'
++       !%filscr_it1=TRIM(filnam_ds(3))//'_SCR'
+        filscr_it1=Dtfil%fnameabi_scr
+      end if
+      !    
+@@ -193,7 +193,7 @@
+      if (read_sus1/=0) then ! Change the default.
+        Dtfil%fnameabi_sus=filchi0_it1
+      else                   ! Save the IT1 name for GW0 runs
+-       !$filchi0_it1=TRIM(filnam_ds(3))//'_SUS'
++       !%filchi0_it1=TRIM(filnam_ds(3))//'_SUS'
+        filchi0_it1 =Dtfil%fnameabi_sus
+      end if
+    end if
+diff -Naur abinit-6.4.2.bak/src/95_drive/iofn1.F90 abinit-6.4.2/src/95_drive/iofn1.F90
+--- src/95_drive/iofn1.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/95_drive/iofn1.F90	2011-01-05 09:17:55.000000000 +0000
+@@ -197,10 +197,10 @@
+  call filnam_comm(5,filnam,mpi_enreg)
+ 
+ !MG TODO
+-!$call xbarrier_mpi(comm_world)
+-!$call timab(48,1,tsec)
+-!$call xcast_mpi(filnam,0,comm_world,ierr)
+-!$call timab(48,2,tsec)
++!%call xbarrier_mpi(comm_world)
++!%call timab(48,1,tsec)
++!%call xcast_mpi(filnam,0,comm_world,ierr)
++!%call timab(48,2,tsec)
+ 
+ !Create a name for the status file, based on filnam(5)
+  filstat=trim(filnam(5))//'_STATUS'
+diff -Naur abinit-6.4.2.bak/src/95_drive/sigma.F90 abinit-6.4.2/src/95_drive/sigma.F90
+--- src/95_drive/sigma.F90	2010-12-01 19:47:18.000000000 +0000
++++ src/95_drive/sigma.F90	2011-01-05 09:17:02.000000000 +0000
+@@ -1505,12 +1505,12 @@
+ !pass it to the setup_ppmodel
+ !* It would be possible to calculate rho(G) using Paw_pwff, though. Maybe faster but
+ !results will depend on the expression used for the matrix elements. This approach is safer.
+-!$ if (Psps%usepaw==1.and.Ppm%needs_rhog>0) then
+-!$  allocate(ks_rhor_paw(nfftf,Dtset%nspden))
+-!$  call denfgr(Cryst%natom,Dtset%nspden,ks_nhat,Cryst%ntypat,Pawfgr,Pawfgrtab,Pawrad,KS_pawrhoij,Pawtab,&
+-!$&   ks_rhor,ks_rhor_paw,Psps,Cryst%typat)
+-!$  deallocate(ks_rhor_paw)
+-!$ end if
++!% if (Psps%usepaw==1.and.Ppm%needs_rhog>0) then
++!%  allocate(ks_rhor_paw(nfftf,Dtset%nspden))
++!%  call denfgr(Cryst%natom,Dtset%nspden,ks_nhat,Cryst%ntypat,Pawfgr,Pawfgrtab,Pawrad,KS_pawrhoij,Pawtab,&
++!%&   ks_rhor,ks_rhor_paw,Psps,Cryst%typat)
++!%  deallocate(ks_rhor_paw)
++!% end if
+ 
+  call timab(409,2,tsec) ! getW
+ !