From 2798183c6d704e78b451cb72e00503d71e448a85 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Sat, 19 Nov 2022 16:16:04 +0100
Subject: [PATCH] Removed labels in loc.f

---
 src/ao_one_e_ints/pseudopot.f90 |   8 +-
 src/utils/loc.f                 | 297 +++++++++++++++++---------------
 2 files changed, 163 insertions(+), 142 deletions(-)

diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90
index 141292d8..bad641ab 100644
--- a/src/ao_one_e_ints/pseudopot.f90
+++ b/src/ao_one_e_ints/pseudopot.f90
@@ -1095,9 +1095,9 @@ double precision function overlap_orb_ylm_grid(nptsgrid,r_orb,npower_orb,center_
 implicit none
 !! PSEUDOS
 integer nptsgridmax,nptsgrid
-double precision coefs_pseudo,ptsgrid
 parameter(nptsgridmax=50)
-common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3)
+double precision coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3)
+common/pseudos/coefs_pseudo,ptsgrid
 !!!!!
 integer npower_orb(3),l,m,i
 double precision x,g_orb,two_pi,dx,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3)
@@ -1235,10 +1235,10 @@ end
       subroutine initpseudos(nptsgrid)
       implicit none
       integer nptsgridmax,nptsgrid,ik
-      double precision coefs_pseudo,ptsgrid
       double precision p,q,r,s
       parameter(nptsgridmax=50)
-      common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3)
+      double precision :: coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3)
+      common/pseudos/coefs_pseudo,ptsgrid
 
       p=1.d0/dsqrt(2.d0)
       q=1.d0/dsqrt(3.d0)
diff --git a/src/utils/loc.f b/src/utils/loc.f
index 6e5d7345..02693281 100644
--- a/src/utils/loc.f
+++ b/src/utils/loc.f
@@ -1,6 +1,6 @@
 c************************************************************************
       subroutine maxovl(n,m,s,t,w)
-C 
+C
 C     This subprogram contains an iterative procedure to find the
 C     unitary transformation of a set of n vectors which maximizes
 C     the sum of their square overlaps with a set of m reference
@@ -10,7 +10,7 @@ C     S: overlap matrix <ref|vec>
 C     T: rotation matrix
 C     W: new overlap matrix
 C
-C 
+C
       implicit real*8(a-h,o-y),logical*1(z)
 !      parameter (id1=700)
 !      dimension s(id1,id1),t(id1,id1),w(id1,id1)
@@ -29,23 +29,26 @@ C      conv=1.d-6
      * 5x,'following the principle of maximum overlap with a set of',
      * i3,' reference vectors'/5x,'required convergence on rotation ',
      * 'angle =',f13.10///5x,'Starting overlap matrix'/)
-      do 6 i=1,m
-      write (6,145) i
-    6 write (6,150) (s(i,j),j=1,n)
+      do i=1,m
+        write (6,145) i
+        write (6,150) (s(i,j),j=1,n)
+      end do
     8 mm=m-1
       if (m.lt.n) mm=m
       iter=0
-      do 20 j=1,n
-      do 16 i=1,n
-      t(i,j)=0.d0
-   16 continue
-      do 18 i=1,m
-   18 w(i,j)=s(i,j)
-   20 t(j,j)=1.d0
+      do j=1,n
+        do i=1,n
+          t(i,j)=0.d0
+        end do
+        do i=1,m
+          w(i,j)=s(i,j)
+        enddo
+        t(j,j)=1.d0
+      enddo
       sum=0.d0
-      do 10 i=1,m
-      sum=sum+s(i,i)*s(i,i)
-   10 continue
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      end do
       sum=sum/m
       if (zprt) write (6,12) sum
    12 format (//5x,'Average square overlap =',f10.6)
@@ -54,18 +57,18 @@ C      conv=1.d-6
       j=1
    21 if (j.ge.last) goto 30
       sum=0.d0
-                                                                
-      do 22 i=1,n
-   22 sum=sum+s(i,j)*s(i,j)
+      do i=1,n
+        sum=sum+s(i,j)*s(i,j)
+      enddo
       if (sum.gt.small) goto 28
-      do 24 i=1,n
-      sij=s(i,j)
-      s(i,j)=-s(i,last)
-      s(i,last)=sij
-      tij=t(i,j)
-      t(i,j)=-t(i,last)
-      t(i,last)=tij
-   24 continue
+      do i=1,n
+        sij=s(i,j)
+        s(i,j)=-s(i,last)
+        s(i,last)=sij
+        tij=t(i,j)
+        t(i,j)=-t(i,last)
+        t(i,last)=tij
+      end do
       last=last-1
       goto 21
    28 j=j+1
@@ -101,17 +104,18 @@ C      conv=1.d-6
       sine=1.d0
    34 delta=sine*(a*sine+b*cosine)
       if (zprt.and.delta.lt.0.d0) write (6,71) i,j,a,b,sine,cosine,delta
-      do 35 k=1,m
-      p=s(k,i)*cosine-s(k,j)*sine
-      q=s(k,i)*sine+s(k,j)*cosine
-      s(k,i)=p
-   35 s(k,j)=q
-      do 40 k=1,n
-      p=t(k,i)*cosine-t(k,j)*sine
-      q=t(k,i)*sine+t(k,j)*cosine
-      t(k,i)=p
-      t(k,j)=q
-   40 continue
+      do k=1,m
+        p=s(k,i)*cosine-s(k,j)*sine
+        q=s(k,i)*sine+s(k,j)*cosine
+        s(k,i)=p
+        s(k,j)=q
+      enddo
+      do k=1,n
+        p=t(k,i)*cosine-t(k,j)*sine
+        q=t(k,i)*sine+t(k,j)*cosine
+        t(k,i)=p
+        t(k,j)=q
+      enddo
    45 d=dabs(sine)
       if (d.le.amax) goto 50
       imax=i
@@ -132,43 +136,50 @@ C      conv=1.d-6
      * 'in subroutine maxovl ***'//)
       stop
   100 continue
-      do 120 j=1,n
-      if (s(j,j).gt.0.d0) goto 120
-      do 105 i=1,m
-  105 s(i,j)=-s(i,j)
-      do 110 i=1,n
-  110 t(i,j)=-t(i,j)
-  120 continue
+      do j=1,n
+        if (s(j,j).gt.0.d0) cycle
+        do i=1,m
+          s(i,j)=-s(i,j)
+        enddo
+        do i=1,n
+          t(i,j)=-t(i,j)
+        enddo
+      enddo
       sum=0.d0
-      do 125 i=1,m
-  125 sum=sum+s(i,i)*s(i,i)
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      enddo
       sum=sum/m
-      do 122 i=1,m
-      do 122 j=1,n
-      sw=s(i,j)
-      s(i,j)=w(i,j)
-  122 w(i,j)=sw
+      do i=1,m
+        do j=1,n
+          sw=s(i,j)
+          s(i,j)=w(i,j)
+          w(i,j)=sw
+        enddo
+      enddo
       if (.not.zprt) return
       write (6,12) sum
       write (6,130)
   130 format (//5x,'transformation matrix')
-      do 140 i=1,n
-      write (6,145) i
-  140 write (6,150) (t(i,j),j=1,n)
+      do i=1,n
+        write (6,145) i
+        write (6,150) (t(i,j),j=1,n)
+      enddo
   145 format (i8)
   150 format (2x,10f12.8)
       write (6,160)
   160 format (//5x,'new overlap matrix'/)
-      do 170 i=1,m
-      write (6,145) i
-  170 write (6,150) (w(i,j),j=1,n)
+      do i=1,m
+        write (6,145) i
+        write (6,150) (w(i,j),j=1,n)
+      enddo
       return
       end
 
 
 c************************************************************************
       subroutine maxovl_no_print(n,m,s,t,w)
-C 
+C
 C     This subprogram contains an iterative procedure to find the
 C     unitary transformation of a set of n vectors which maximizes
 C     the sum of their square overlaps with a set of m reference
@@ -178,7 +189,7 @@ C     S: overlap matrix <ref|vec>
 C     T: rotation matrix
 C     W: new overlap matrix
 C
-C 
+C
       implicit real*8(a-h,o-y),logical*1(z)
       parameter (id1=300)
       dimension s(id1,id1),t(id1,id1),w(id1,id1)
@@ -193,17 +204,19 @@ C      conv=1.d-6
     8 mm=m-1
       if (m.lt.n) mm=m
       iter=0
-      do 20 j=1,n
-      do 16 i=1,n
-      t(i,j)=0.d0
-   16 continue
-      do 18 i=1,m
-   18 w(i,j)=s(i,j)
-   20 t(j,j)=1.d0
+      do j=1,n
+        do i=1,n
+          t(i,j)=0.d0
+        enddo
+        do i=1,m
+          w(i,j)=s(i,j)
+        enddo
+        t(j,j)=1.d0
+      enddo
       sum=0.d0
-      do 10 i=1,m
-      sum=sum+s(i,i)*s(i,i)
-   10 continue
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      enddo
       sum=sum/m
    12 format (//5x,'Average square overlap =',f10.6)
       if (n.eq.1) goto 100
@@ -211,18 +224,19 @@ C      conv=1.d-6
       j=1
    21 if (j.ge.last) goto 30
       sum=0.d0
-                                                                
-      do 22 i=1,n
-   22 sum=sum+s(i,j)*s(i,j)
+
+      do i=1,n
+        sum=sum+s(i,j)*s(i,j)
+      enddo
       if (sum.gt.small) goto 28
-      do 24 i=1,n
-      sij=s(i,j)
-      s(i,j)=-s(i,last)
-      s(i,last)=sij
-      tij=t(i,j)
-      t(i,j)=-t(i,last)
-      t(i,last)=tij
-   24 continue
+      do i=1,n
+        sij=s(i,j)
+        s(i,j)=-s(i,last)
+        s(i,last)=sij
+        tij=t(i,j)
+        t(i,j)=-t(i,last)
+        t(i,last)=tij
+      end do
       last=last-1
       goto 21
    28 j=j+1
@@ -232,50 +246,52 @@ C      conv=1.d-6
       jmax=0
       dmax=0.d0
       amax=0.d0
-      do 60 i=1,mm
-      ip=i+1
-      do 50 j=ip,n
-      a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
-      b=-s(i,i)*s(i,j)
-      if (j.gt.m) goto 31
-      a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
-      b=b+s(j,i)*s(j,j)
-   31 b=b+b
-      if (a.eq.0.d0) goto 32
-      ba=b/a
-      if (dabs(ba).gt.small) goto 32
-      if (a.gt.0.d0) goto 33
-      tang=-0.5d0*ba
-      cosine=1.d0/dsqrt(1.d0+tang*tang)
-      sine=tang*cosine
-      goto 34
-   32 tang=0.d0
-      if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
-      cosine=1.d0/dsqrt(1.d0+tang*tang)
-      sine=tang*cosine
-      goto 34
-   33 cosine=0.d0
-      sine=1.d0
-   34 delta=sine*(a*sine+b*cosine)
-      do 35 k=1,m
-      p=s(k,i)*cosine-s(k,j)*sine
-      q=s(k,i)*sine+s(k,j)*cosine
-      s(k,i)=p
-   35 s(k,j)=q
-      do 40 k=1,n
-      p=t(k,i)*cosine-t(k,j)*sine
-      q=t(k,i)*sine+t(k,j)*cosine
-      t(k,i)=p
-      t(k,j)=q
-   40 continue
-   45 d=dabs(sine)
-      if (d.le.amax) goto 50
-      imax=i
-      jmax=j
-      amax=d
-      dmax=delta
-   50 continue
-   60 continue
+      do i=1,mm
+        ip=i+1
+        do j=ip,n
+          a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
+          b=-s(i,i)*s(i,j)
+          if (j.gt.m) goto 31
+          a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
+          b=b+s(j,i)*s(j,j)
+  31      b=b+b
+          if (a.eq.0.d0) goto 32
+          ba=b/a
+          if (dabs(ba).gt.small) goto 32
+          if (a.gt.0.d0) goto 33
+          tang=-0.5d0*ba
+          cosine=1.d0/dsqrt(1.d0+tang*tang)
+          sine=tang*cosine
+          goto 34
+  32      tang=0.d0
+          if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
+          cosine=1.d0/dsqrt(1.d0+tang*tang)
+          sine=tang*cosine
+          goto 34
+  33      cosine=0.d0
+          sine=1.d0
+  34      delta=sine*(a*sine+b*cosine)
+          do k=1,m
+            p=s(k,i)*cosine-s(k,j)*sine
+            q=s(k,i)*sine+s(k,j)*cosine
+            s(k,i)=p
+            s(k,j)=q
+          enddo
+          do k=1,n
+            p=t(k,i)*cosine-t(k,j)*sine
+            q=t(k,i)*sine+t(k,j)*cosine
+            t(k,i)=p
+            t(k,j)=q
+          enddo
+  45      d=dabs(sine)
+          if (d.le.amax) goto 50
+          imax=i
+          jmax=j
+          amax=d
+          dmax=delta
+  50      continue
+        end do
+      end do
    70 format (' iter=',i4,' largest rotation=',f12.8,
      * ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5)
    71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5)
@@ -285,22 +301,27 @@ C      conv=1.d-6
      * 'in subroutine maxovl ***'//)
       stop
   100 continue
-      do 120 j=1,n
-      if (s(j,j).gt.0.d0) goto 120
-      do 105 i=1,m
-  105 s(i,j)=-s(i,j)
-      do 110 i=1,n
-  110 t(i,j)=-t(i,j)
-  120 continue
+      do j=1,n
+        if (s(j,j).gt.0.d0) cycle
+        do i=1,m
+          s(i,j)=-s(i,j)
+        enddo
+        do i=1,n
+          t(i,j)=-t(i,j)
+        enddo
+      enddo
       sum=0.d0
-      do 125 i=1,m
-  125 sum=sum+s(i,i)*s(i,i)
+      do i=1,m
+        sum=sum+s(i,i)*s(i,i)
+      enddo
       sum=sum/m
-      do 122 i=1,m
-      do 122 j=1,n
-      sw=s(i,j)
-      s(i,j)=w(i,j)
-  122 w(i,j)=sw
+      do i=1,m
+        do j=1,n
+          sw=s(i,j)
+          s(i,j)=w(i,j)
+          w(i,j)=sw
+        enddo
+      enddo
       return
       end