Remove workarounds

src/dft.f90

@@ -36,8 +36,8 @@ subroutine V2rho(d)
 ! Calculates rho from V by solving Kohn-Sham equations
 type(dft_data_t), intent(inout) :: d
 
-real(dp), dimension(size(d%R)) :: P, Q, Y, tmp
-integer :: converged, i, n, l, relat, j
+real(dp), dimension(size(d%R)) :: P, Q, Y
+integer :: converged, i, n, l, relat
 real(dp) :: Ein, Emin_init, Emax_init
 
 d%rho(:) = 0
@@ -75,9 +75,7 @@ subroutine V2rho(d)
     else
         Y = sqrt(P**2 + Q**2) / d%R
     end if
-    do j = 1, size(d%R)
-        d%rho(j) = d%rho(j) + d%fo(i) * Y(j)**2
-    end do
+    d%rho = d%rho + d%fo(i) * Y**2
     d%orbitals(:, i) = Y
 end do
 d%rho = d%rho / (4*pi)
@@ -109,31 +107,19 @@ subroutine total_energy(fo, ks_energies, V_in, V_h, V_coulomb, e_xc, R, Rp, n, &
 real(dp), intent(in) :: R(:), Rp(:), n(:) ! Radial grid, number density (positive)
 real(dp), intent(out) :: Etot ! Total energy
 real(dp), intent(out) :: T_s, E_ee, E_en, EE_xc ! Parts of the total energy
-real(dp) :: tmp(size(R))
+
 real(dp) :: rho(size(n))
-integer :: i
 real(dp) :: E_c, E_band!, Exc2
 rho = -n
 
 E_band = sum(fo * ks_energies)
-do i = 1, size(R)
-    tmp(i) = V_in(i) * rho(i) * R(i)**2
-end do
-T_s = E_band + 4*pi * integrate(Rp, tmp)
+T_s = E_band + 4*pi * integrate(Rp, V_in * rho * R**2)
 
-do i = 1, size(R)
-    tmp(i) = V_h(i) * rho(i) * R(i)**2
-end do
-E_ee = -2*pi * integrate(Rp, tmp)
-do i = 1, size(R)
-    tmp(i) = (-V_coulomb(i)) * rho(i) * R(i)**2
-end do
-E_en =  4*pi * integrate(Rp, tmp)
+E_ee = -2*pi * integrate(Rp, V_h * rho * R**2)
+E_en =  4*pi * integrate(Rp, (-V_coulomb) * rho * R**2)
 E_c = E_ee + E_en
-do i = 1, size(R)
-    tmp(i) = e_xc(i) * rho(i) * R(i)**2
-end do
-EE_xc = -4*pi * integrate(Rp, tmp)
+
+EE_xc = -4*pi * integrate(Rp, e_xc * rho * R**2)
 
 Etot = T_s + E_c + EE_xc
 end subroutine

src/energies.f90

@@ -48,7 +48,6 @@ function thomas_fermi_potential(R, Z, cut) result(V)
 logical, intent(in), optional :: cut ! Cut the potential, default .true.
 real(dp) :: x(size(R)), Z_eff(size(R)), V(size(R))
 real(dp) :: alpha, beta, gamma
-integer :: i
 
 x = R * (128*Z/(9*pi**2)) ** (1.0_dp/3)
 ! Z_eff(x) = Z * phi(x), where phi(x) satisfies the Thomas-Fermi equation:
@@ -64,11 +63,7 @@ function thomas_fermi_potential(R, Z, cut) result(V)
 Z_eff = Z * (1 + alpha*sqrt(x) + beta*x*exp(-gamma*sqrt(x)))**2 * &
     exp(-2*alpha*sqrt(x))
 ! This keeps all the eigenvalues of the radial problem negative:
-if (.not. present(cut)) then
-    do i = 1, size(R)
-        if( Z_eff(i) < 1 ) Z_eff(i) = 1
-    end do
-end if
+if (.not. present(cut)) where (Z_eff < 1) Z_eff = 1
 V = -Z_eff / r
 end function
 

src/ode1d.f90

@@ -163,20 +163,14 @@ subroutine parsefunction(y, nodes, minidx, positive)
 integer function get_n_nodes(y) result(nodes)
 ! Returns the number of nodes of the function 'y'
 real(dp), intent(in) :: y(:)
-integer :: sizey
 
 integer :: last_sign, last_i, i, isy
 nodes = 0
 last_sign = int(sign(1.0_dp, y(1)))
 last_i = -1
-sizey = size(y)
 
-do i = 2, sizey
-  if( y(i) >= 0 ) then
-    isy = 1
-  else
-    isy = -1
-  end if
+do i = 2, size(y)
+  isy = int(sign(1.0_dp, y(i)))
   if (isy == -last_sign) then
       last_sign = isy
       last_i = i - 1

src/rdirac.f90

@@ -45,27 +45,21 @@ subroutine dirac_outward_adams(c, kappa, Z, E, R, Rp, V, P, Q, imax)
 integer :: i
 real(dp) :: lambda, Delta, M(2, 2), u1_tmp, u2_tmp
 real(dp) :: Vmid(3)
-real(dp) :: tmp(4)
 
 if (size(R) < 4) call stop_error("size(R) <= 4")
-tmp = R(:4)
-Vmid = get_midpoints(tmp, V(:4))
-! Vmid = get_midpoints(R(:4), V(:4))
+Vmid = get_midpoints(R(:4), V(:4))
 call integrate_radial_dirac_r_rk4(c, kappa, Z, E, R(:4), V(:4), Vmid, &
     u1(:4), u2(:4), imax)
 if (imax /= 4) call stop_error("rk4 failed")
 
-do i = 1, size(R)
-    Ctot(i, 1, 1) = -kappa / R(i)
-    Ctot(i, 2, 2) = +kappa / R(i)
-    Ctot(i, 1, 2) = +(E-V(i))/c + 2*c
-    Ctot(i, 2, 1) = -(E-V(i))/c
-end do
+Ctot(:, 1, 1) = -kappa / R
+Ctot(:, 2, 2) = +kappa / R
+Ctot(:, 1, 2) = +(E-V)/c + 2*c
+Ctot(:, 2, 1) = -(E-V)/c
+
+u1p(:4) = Rp(:4) * (Ctot(:4, 1, 1) * u1(:4) + Ctot(:4, 1, 2) * u2(:4))
+u2p(:4) = Rp(:4) * (Ctot(:4, 2, 1) * u1(:4) + Ctot(:4, 2, 2) * u2(:4))
 
-do i = 1, 4
-    u1p(i) = Rp(i) * (Ctot(i, 1, 1) * u1(i) + Ctot(i, 1, 2) * u2(i))
-    u2p(i) = Rp(i) * (Ctot(i, 2, 1) * u1(i) + Ctot(i, 2, 2) * u2(i))
-end do
 do i = 4, size(R)-1
     u1p(i) = Rp(i) * (Ctot(i, 1, 1)*u1(i) + Ctot(i, 1, 2)*u2(i))
     u2p(i) = Rp(i) * (Ctot(i, 2, 1)*u1(i) + Ctot(i, 2, 2)*u2(i))
@@ -149,24 +143,20 @@ subroutine dirac_inward_adams(c, kappa, E, R, Rp, V, P, Q, imin)
 u1(i_max+4:) = 0
 u2(i_max+4:) = 0
 
-do i = i_max, i_max+4
-    u1(i) = exp(-lambda * (R(i)-R(1))) / sqrt(-E/(E+2*c**2))
-    u2(i) = -exp(-lambda * (R(i)-R(1)))
-end do
-do i = 1, size(R)
-    Ctot(i, 1, 1) = -kappa / R(i)
-    Ctot(i, 2, 2) = +kappa / R(i)
-    Ctot(i, 1, 2) = +(E-V(i))/c + 2*c
-    Ctot(i, 2, 1) = -(E-V(i))/c
-end do
-do i = i_max, i_max+4
-    u1p(i) = Rp(i) * &
-        (Ctot(i, 1, 1)*u1(i) &
-            + Ctot(i, 1, 2) * u2(i))
-    u2p(i) = Rp(i) * &
-        (Ctot(i, 2, 1)*u1(i) &
-            + Ctot(i, 2, 2) * u2(i))
-end do
+u1(i_max:i_max+4) = exp(-lambda * (R(i_max:i_max+4)-R(1))) / sqrt(-E/(E+2*c**2))
+u2(i_max:i_max+4) = -exp(-lambda * (R(i_max:i_max+4)-R(1)))
+
+Ctot(:, 1, 1) = -kappa / R
+Ctot(:, 2, 2) = +kappa / R
+Ctot(:, 1, 2) = +(E-V)/c + 2*c
+Ctot(:, 2, 1) = -(E-V)/c
+
+u1p(i_max:i_max+4) = Rp(i_max:i_max+4) * &
+    (Ctot(i_max:i_max+4, 1, 1)*u1(i_max:i_max+4) &
+        + Ctot(i_max:i_max+4, 1, 2) * u2(i_max:i_max+4))
+u2p(i_max:i_max+4) = Rp(i_max:i_max+4) * &
+    (Ctot(i_max:i_max+4, 2, 1)*u1(i_max:i_max+4) &
+        + Ctot(i_max:i_max+4, 2, 2) * u2(i_max:i_max+4))
 
 do i = i_max, 2, -1
     u1p(i) = Rp(i) * (Ctot(i, 1, 1)*u1(i) + Ctot(i, 1, 2)*u2(i))
@@ -292,7 +282,7 @@ subroutine integrate_radial_dirac_r_rk4(c, kappa, Z, E, R, V, Vmid, P, Q, imax)
 real(dp), dimension(size(R)-1) :: Rmid
 
 real(dp) :: beta, Z1
-integer :: l, i
+integer :: l
 
 beta = sqrt(kappa**2-(Z/c)**2)
 if (Z /= 0) then
@@ -311,21 +301,16 @@ subroutine integrate_radial_dirac_r_rk4(c, kappa, Z, E, R, V, Vmid, P, Q, imax)
     end if
 end if
 
-do i = 1, size(R)
-    Ctot(i, 1, 1) = -kappa / R(i)
-    Ctot(i, 2, 2) = +kappa / R(i)
-    Ctot(i, 1, 2) = +(E-V(i))/c + 2*c
-    Ctot(i, 2, 1) = -(E-V(i))/c
-end do
-do i = 1, size(Rmid)
-    Rmid(i) = (R(i) + R(i+1)) / 2
-end do
-do i = 1, size(R) - 1
-    Cmid(i, 1, 1) = -kappa / Rmid(i)
-    Cmid(i, 2, 2) = +kappa / Rmid(i)
-    Cmid(i, 1, 2) = +(E-Vmid(i))/c + 2*c
-    Cmid(i, 2, 1) = -(E-Vmid(i))/c
-end do
+Ctot(:, 1, 1) = -kappa / R
+Ctot(:, 2, 2) = +kappa / R
+Ctot(:, 1, 2) = +(E-V)/c + 2*c
+Ctot(:, 2, 1) = -(E-V)/c
+
+Rmid = (R(:size(R)-1) + R(2:)) / 2
+Cmid(:, 1, 1) = -kappa / Rmid
+Cmid(:, 2, 2) = +kappa / Rmid
+Cmid(:, 1, 2) = +(E-Vmid)/c + 2*c
+Cmid(:, 2, 1) = -(E-Vmid)/c
 call rk4_integrate4(R, y, Ctot, Cmid, max_val, P, Q, imax)
 end subroutine
 

src/reigen.f90

@@ -169,13 +169,15 @@ subroutine solve_radial_eigenproblem(n, l, Ein, eps, max_iter, &
 E = Ein
 if (.not.(n > 0)) call stop_error("n > 0 not satisfied")
 if (.not.((0 <= l).and.(l < n))) call stop_error("0 <= l < n not satisfied")
+
 Emax = Emax_init
 Emin = Emin_init
 if (E > Emax .or. E < Emin) E = (Emin + Emax) / 2
 
 iter = 0
 last_bisect = .true.
-l1: do iter = 0, max_iter - 1
+do while (iter < max_iter)
+    iter = iter + 1
 
     ! See if bisection is converged
     if (abs(Emax - Emin) < eps) then
@@ -258,7 +260,7 @@ subroutine solve_radial_eigenproblem(n, l, Ein, eps, max_iter, &
         end if
         E = (Emin + Emax) / 2
         last_bisect = .true.
-        cycle l1
+        cycle
     end if
 
     ! Perturbation theory correction
@@ -312,7 +314,7 @@ subroutine solve_radial_eigenproblem(n, l, Ein, eps, max_iter, &
         E = E + dE
         last_bisect = .false.
     end if
-end do l1
+end do
 if (iter == max_iter) then
     ! We didn't converge in 'max_iter' iterations
     converged = 2

src/rpoisson.f90

@@ -37,12 +37,9 @@ function rpoisson_outward_pc(R, Rp, rho) result(V)
 integer :: N, i, it
 integer, parameter :: max_it = 2
 real(dp) :: rho_mid(3)
-real(dp) :: tmp(4)
 
 N = size(R)
-tmp = R(:4)
-rho_mid = get_midpoints(tmp, rho(:4))
-! rho_mid = get_midpoints(R(:4), rho(:4))
+rho_mid = get_midpoints(R(:4), rho(:4))
 call rpoisson_outward_rk4(rho(:4), rho_mid, R(:4), &
     4*pi*integrate(Rp, rho*R), &
     0.0_dp, &

src/rschroed.f90

@@ -54,17 +54,14 @@ subroutine schroed_outward_adams(l, Z, E, R, Rp, V, P, Q, imax)
 real(dp), dimension(size(R)) :: C, u1, u2, u1p, u2p, Vmid
 integer :: i
 real(dp) :: lambda, Delta, M(2, 2), u1_tmp, u2_tmp
-real(dp) :: tmp(4)
 
 if (size(R) < 5) call stop_error("size(R) < 5")
 if (.not. (size(R) == size(Rp) .and. size(R) == size(V) .and. &
     size(R) == size(P) .and. size(R) == size(P) .and. size(R) == size(Q))) then
     call stop_error("Array sizes mismatch")
 end if
 C = (l*(l+1)/R**2 + 2 * (V-E))
-tmp = R(:4)
-Vmid(:3) = get_midpoints(tmp, V(:4))
-! Vmid(:3) = get_midpoints(R(:4), V(:4))
+Vmid(:3) = get_midpoints(R(:4), V(:4))
 call integrate_rschroed_rk4(l, Z, E, R(:4), V(:4), Vmid(:3), &
     u1(:4), u2(:4), imax)
 !u1(1:4) = R(1:4) ** (l+1)
@@ -124,7 +121,6 @@ subroutine schroed_inward_adams(l, E, R, Rp, V, P, Q, imin)
 integer :: i, i_max
 real(dp) :: lambda, Delta, M(2, 2), u1_tmp, u2_tmp
 real(dp) :: R_max
-integer :: imaxitr
 
 C = (l*(l+1)/R**2 + 2 * (V-E))
 
@@ -144,12 +140,10 @@ subroutine schroed_inward_adams(l, E, R, Rp, V, P, Q, imin)
     i_max = i_max - 1
 end do
 
-do imaxitr = i_max, i_max + 4
-    u1(imaxitr) = exp(-lambda * (R(imaxitr)-R(1)))
-    u2(imaxitr) = - lambda * u1(imaxitr)
-    u1p(imaxitr) = Rp(imaxitr) * u2(imaxitr)
-    u2p(imaxitr) = Rp(imaxitr) * C(imaxitr) * u1(imaxitr)
-end do
+u1(i_max:i_max+4) = exp(-lambda * (R(i_max:i_max+4)-R(1)))
+u2(i_max:i_max+4) = - lambda * u1(i_max:i_max+4)
+u1p(i_max:i_max+4) = Rp(i_max:i_max+4)                * u2(i_max:i_max+4)
+u2p(i_max:i_max+4) = Rp(i_max:i_max+4) * C(i_max:i_max+4) * u1(i_max:i_max+4)
 
 do i = i_max, 2, -1
     u1p(i) = Rp(i)        * u2(i)