Slight cleanup of bfgs*

force/CMakeLists.txt

 set(fleur_F77 ${fleur_F77}
-force/bfgs.f
-force/bfgs0.f
 force/relax.F
 )
 set(fleur_F90 ${fleur_F90}
+force/bfgs0.f90
+force/bfgs.f90
 force/force_b8.f90
 force/force_0.f90
 force/force_a12.f90

force/bfgs.f

-      MODULE m_bfgs
-      use m_juDFT
-      CONTAINS
-      SUBROUTINE bfgs(
-     >                ntype,istep,istep0,force,
-     >                zat,xa,thetad,epsdisp,epsforce,tote,
-     X                xold,y,h,tau0,
-     <                lconv)
-!*******************************************************************
-!     relaxes the forces using the BFGS quasi-Newton method.
-!     input:
-!     istep  = atomic step in this run
-!     istep0 = number of atomic steps in previous runs used
-!              in updating hessian
-!     
-!     output:
-!        lconv   = logical true if forces are converged to tolerance
-!                  given in epsforce
-!     
-!     the positions and forces from this step are added 
-!     to file force.dat
-!*******************************************************************
-      
-      IMPLICIT NONE
-!     ..
-!     .. Scalar Arguments ..
-      INTEGER, INTENT (IN)    :: ntype,istep,istep0
-      REAL,    INTENT (IN)    :: epsdisp,epsforce,tote
-      REAL,    INTENT (INOUT) :: thetad,xa
-      LOGICAL, INTENT (OUT)   :: lconv
-!     ..
-!     .. Array Arguments ..
-      REAL,    INTENT (IN)    :: force(3,ntype),zat(ntype)
-      REAL,    INTENT (INOUT) :: tau0(3,ntype)
-      REAL,    INTENT (INOUT) :: y(3*ntype),xold(3*ntype)
-      REAL,    INTENT (INOUT) :: h(3*ntype,3*ntype)
-!     ..
-!     .. Local Scalars ..
-      INTEGER :: i,j,n,nn,ist,na
-      REAL    :: py,yy,alpha,s,zatm,fmax,gamma,d2,dispmax
-!     ..
-!     .. Local  Arrays ..
-      REAL :: f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype)
-!     
-      n=3*ntype
-      ist = istep+istep0
-      
-!---  >       get positions and forces in correct form and output
-      DO na = 1,ntype
-         nn = 3*(na-1)
-         xnew(nn+1:nn+3) = tau0(:,na)
-         f(nn+1:nn+3)    = force(:,na)
-      ENDDO
-      !Write new entry into forces.dat
-      OPEN(43,file ='forces.dat',status ='unknown',form='formatted'
-     $     ,position ='append')
-      WRITE (43,'(a,f20.10)') "energy =",tote
-      WRITE (43,'(3f16.9,3f14.9)') 
-     >     ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1)
-      CLOSE(43)
-!---  >       get maximum force
-      fmax = 0.0
-      DO na = 1,ntype
-         nn = 3*(na-1)
-         fmax = MAX( fmax, (f(nn+1)**2+f(nn+2)**2+f(nn+3)**2) )
-      ENDDO
-      fmax = SQRT(fmax)
-      WRITE (6,1000) istep,fmax
-      IF ( fmax<epsforce) THEN
-         lconv = .TRUE.
-         RETURN
-      ELSE
-         lconv = .FALSE.
-      ENDIF
- 1000 FORMAT (1x/,' atomic step',i4,': maximum force =',
-     +     1p,e14.6,' hartrees/a.u.')
-
-!------------------------------------------------------------
-! ===>       if first step, go along gradient
-      
-      IF (ist==1) THEN
-!---  >       choose a reasonable first guess for scaling, but
-!---  >       limit displacement to a maximum of 0.25 a.u.
-!---  >       (may need to be changed for different systems)
-!---  >       this choice is based on a Debye temperature of 330K;
-!---  >       modify as needed (change thetad in param.8)
-         zatm = 0.0
-         DO i = 1,ntype
-            zatm = MAX(zatm,zat(i))
-         ENDDO
-         IF (ABS(xa)<1.0e-10) THEN
-            WRITE (6,*) 'WARNING, xa = 0.0 set to 2.0'
-            xa = 2.0
-         ENDIF
-         IF (ABS(thetad)<1.0e-10) THEN
-            WRITE (6,*) 'WARNING, thetad = 0.0 set to 330.0'
-            thetad = 330.0
-         ENDIF
-         alpha = (250.0/(zatm*xa))*((330./thetad)**2)
-         IF ( alpha*fmax*xa > 0.15 ) alpha = 0.25/(fmax*xa)
-         p(:) = alpha*f(:)
-!---  >       set h to identity matrix
-         h = 0.0
-         DO j = 1,n
-            h(j,j) = 1.0
-         ENDDO
-         
-      ELSE
-!------------------------------------------------------------
-! ===>              k-th step
-!     
-!---  >       determine p
-         p(:) = xnew(:)-xold(:)
-!---  >       update the change in gradients
-         y(:) = y(:)-f(:)
-!---  >       get necessary inner products and H|y>
-         py = dot_PRODUCT(p,y)
-         yy = 0.0
-         DO i = 1,n
-            s = 0.0
-            DO j = 1,n
-               s = s+y(j)*h(j,i)
-            ENDDO
-            v(i) = s
-            yy   = yy+y(i)*s
-         ENDDO
-!---  >       check that update will leave h positive definite;
-!---  >       if not, then stop
-         IF (py<=0.0) THEN
-            WRITE (6,*) '  bfgs: <p|y> < 0'
-            WRITE (6,*) '  check convergence of forces'
-             CALL juDFT_error("bfgs: <p|y><0",calledby="bfgs")
-         ELSE
-!---  >          update h
-            IF (ist==2) THEN
-               gamma = py/yy
-            ELSE
-               gamma = 1.0
-            ENDIF
-            DO j = 1,n
-               DO i = 1,n
-                  h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma
-     +                 + (1.+gamma*yy/py)*p(i)*p(j)/py
-               ENDDO
-            ENDDO
-!---  >          generate p
-            DO i = 1,n
-               s = 0.0
-               DO j = 1,n
-                  s = s+f(j)*h(j,i)
-               ENDDO
-               p(i) = s
-            ENDDO
-         ENDIF
-      ENDIF
-      
-!-------------------------------------------------------------
-!---  >       put xold and y in the correct form for the next step
-      DO i = 1,n
-         xold(i) = xnew(i)
-         y(i) = f(i)
-      ENDDO
-      
-!---  >    if displacements are all less than epsdisp, then converged
-      dispmax = 0.0
-      DO na = 1,ntype
-         nn = 3*(na-1)
-         d2 = p(nn+1)**2 + p(nn+2)**2 + p(nn+3)**2
-         dispmax = MAX( dispmax, d2)
-      ENDDO
-      dispmax = xa*SQRT(dispmax)
-      IF (dispmax<epsdisp) THEN
-         lconv = .TRUE.
-      ELSE
-         lconv = .FALSE.
-      ENDIF
-      
-!---  >    get new displacements
-      DO i = 1,n
-         xnew(i) = xold(i)+p(i)
-      ENDDO
-      DO na = 1,ntype
-         nn = 3*(na-1)
-         tau0(1,na) = xnew(nn+1)
-         tau0(2,na) = xnew(nn+2)
-         tau0(3,na) = xnew(nn+3)
-      ENDDO
-      
-      WRITE (6,'(1x/)')
-      WRITE (6,*) 'changes in p for step',ist
-      DO na = 1,ntype
-         nn = 3*(na-1)
-         WRITE (6,'(i5,6f12.6)') na,(p(nn+i),i = 1,3),(xnew(nn+j),j=1,3)
-      ENDDO
-      
-      END SUBROUTINE bfgs
-      END MODULE m_bfgs

force/bfgs.f90

+MODULE m_bfgs
+
+USE m_juDFT
+
+CONTAINS
+
+!*******************************************************************
+!     relaxes the forces using the BFGS quasi-Newton method.
+!     input:
+!     istep  = atomic step in this run
+!     istep0 = number of atomic steps in previous runs used
+!              in updating hessian
+!     
+!     output:
+!        lconv   = logical true if forces are converged to tolerance
+!                  given in epsforce
+!     
+!     the positions and forces from this step are added 
+!     to file force.dat
+!*******************************************************************
+SUBROUTINE bfgs(ntype,istep,istep0,force,zat,xa,thetad,epsdisp,epsforce,tote,&
+                xold,y,h,tau0,lconv)
+      
+IMPLICIT NONE
+
+   ! Scalar Arguments
+   INTEGER, INTENT (IN)    :: ntype,istep,istep0
+   REAL,    INTENT (IN)    :: epsdisp,epsforce,tote
+   REAL,    INTENT (INOUT) :: thetad,xa
+   LOGICAL, INTENT (OUT)   :: lconv
+
+   ! Array Arguments
+   REAL,    INTENT (IN)    :: force(3,ntype),zat(ntype)
+   REAL,    INTENT (INOUT) :: tau0(3,ntype)
+   REAL,    INTENT (INOUT) :: y(3*ntype),xold(3*ntype)
+   REAL,    INTENT (INOUT) :: h(3*ntype,3*ntype)
+
+   ! Local Scalars
+   INTEGER :: i,j,n,nn,ist,na
+   REAL    :: py,yy,alpha,s,zatm,fmax,gamma,d2,dispmax
+
+   ! Local Arrays
+   REAL :: f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype)
+
+   n=3*ntype
+   ist = istep+istep0
+      
+   ! get positions and forces in correct form and output
+   DO na = 1,ntype
+      nn = 3*(na-1)
+      xnew(nn+1:nn+3) = tau0(:,na)
+      f(nn+1:nn+3)    = force(:,na)
+   END DO
+
+   ! Write new entry into forces.dat
+   OPEN(43,file ='forces.dat',status ='unknown',form='formatted',position ='append')
+   WRITE (43,'(a,f20.10)') "energy =",tote
+   WRITE (43,'(3f16.9,3f14.9)') ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1)
+   CLOSE(43)
+
+   ! get maximum force
+   fmax = 0.0
+   DO na = 1,ntype
+      nn = 3*(na-1)
+      fmax = MAX( fmax, (f(nn+1)**2+f(nn+2)**2+f(nn+3)**2) )
+   END DO
+   fmax = SQRT(fmax)
+   WRITE (6,1000) istep,fmax
+   IF (fmax < epsforce) THEN
+      lconv = .TRUE.
+      RETURN
+   ELSE
+      lconv = .FALSE.
+   END IF
+   1000 FORMAT (1x/,' atomic step',i4,': maximum force =',1p,e14.6,' hartrees/a.u.')
+
+   ! if first step, go along gradient      
+   IF (ist==1) THEN
+      ! choose a reasonable first guess for scaling, but
+      ! limit displacement to a maximum of 0.25 a.u.
+      ! (may need to be changed for different systems)
+      ! this choice is based on a Debye temperature of 330K;
+      ! modify as needed (change thetad in param.8)
+      zatm = 0.0
+      DO i = 1,ntype
+         zatm = MAX(zatm,zat(i))
+      END DO
+      IF (ABS(xa)<1.0e-10) THEN
+         WRITE (6,*) 'WARNING, xa = 0.0 set to 2.0'
+         xa = 2.0
+      END IF
+      IF (ABS(thetad)<1.0e-10) THEN
+         WRITE (6,*) 'WARNING, thetad = 0.0 set to 330.0'
+         thetad = 330.0
+      END IF
+      alpha = (250.0/(zatm*xa))*((330./thetad)**2)
+      IF ( alpha*fmax*xa > 0.15 ) alpha = 0.25/(fmax*xa)
+      p(:) = alpha*f(:)
+      ! set h to identity matrix
+      h = 0.0
+      DO j = 1,n
+         h(j,j) = 1.0
+      END DO
+   ELSE
+      ! k-th step
+
+      ! determine p
+      p(:) = xnew(:) - xold(:)
+      ! update the change in gradients
+      y(:) = y(:) - f(:)
+      ! get necessary inner products and H|y>
+      py = dot_PRODUCT(p,y)
+      yy = 0.0
+      DO i = 1,n
+         s = 0.0
+         DO j = 1,n
+            s = s+y(j)*h(j,i)
+         END DO
+         v(i) = s
+         yy   = yy+y(i)*s
+      END DO
+      ! check that update will leave h positive definite;
+      ! if not, then stop
+      IF (py <= 0.0) THEN
+         WRITE (6,*) '  bfgs: <p|y> < 0'
+         WRITE (6,*) '  check convergence of forces'
+         CALL juDFT_error("bfgs: <p|y><0",calledby="bfgs")
+      ELSE
+         ! update h
+         IF (ist==2) THEN
+            gamma = py/yy
+         ELSE
+            gamma = 1.0
+         ENDIF
+         DO j = 1,n
+            DO i = 1,n
+               h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma + (1.+gamma*yy/py)*p(i)*p(j)/py
+            END DO
+         END DO
+         ! generate p
+         DO i = 1,n
+            s = 0.0
+            DO j = 1,n
+               s = s+f(j)*h(j,i)
+            END DO
+            p(i) = s
+         END DO
+      END IF
+   END IF
+      
+   ! put xold and y in the correct form for the next step
+   DO i = 1, n
+      xold(i) = xnew(i)
+      y(i) = f(i)
+   END DO
+      
+   ! if displacements are all less than epsdisp, then converged
+   dispmax = 0.0
+   DO na = 1,ntype
+      nn = 3*(na-1)
+      d2 = p(nn+1)**2 + p(nn+2)**2 + p(nn+3)**2
+      dispmax = MAX( dispmax, d2)
+   END DO
+   dispmax = xa*SQRT(dispmax)
+   IF (dispmax<epsdisp) THEN
+      lconv = .TRUE.
+   ELSE
+      lconv = .FALSE.
+   END IF
+
+   ! get new displacements
+   DO i = 1,n
+      xnew(i) = xold(i) + p(i)
+   END DO
+   DO na = 1,ntype
+      nn = 3*(na-1)
+      tau0(1,na) = xnew(nn+1)
+      tau0(2,na) = xnew(nn+2)
+      tau0(3,na) = xnew(nn+3)
+   END DO
+      
+   WRITE (6,'(1x/)')
+   WRITE (6,*) 'changes in p for step',ist
+   DO na = 1,ntype
+      nn = 3*(na-1)
+      WRITE (6,'(i5,6f12.6)') na,(p(nn+i),i = 1,3),(xnew(nn+j),j=1,3)
+   END DO
+
+END SUBROUTINE bfgs
+
+END MODULE m_bfgs

force/bfgs0.f

-      MODULE m_bfgs0
-      CONTAINS
-      SUBROUTINE bfgs0(
-     >                 ntype,
-     <                 istep0,xold,y,h)
-c
-c*******************************************************************
-c     checks whether a file exists with positions and forces to
-c     use to relax structure using the BFGS quasi-Newton method.
-c     istep0 is the number of steps included on file.
-c     h,y, and xold will be saved correctly for further use.
-c*******************************************************************
-c      
-      IMPLICIT NONE
-C     ..
-C     .. Scalar Arguments ..
-      INTEGER, INTENT (IN) :: ntype
-      INTEGER, INTENT (OUT):: istep0
-C     ..
-C     .. Array Arguments ..
-      REAL,    INTENT (OUT):: y(3*ntype),xold(3*ntype)
-      REAL,    INTENT (OUT):: h(3*ntype,3*ntype)
-C     ..
-C     .. Local Scalars ..
-      INTEGER n,is,i,j
-      REAL yy,py,s,gamma
-      CHARACTER :: testchar
-      INTEGER   :: maxhist
-      LOGICAL   :: l_skip
-!     ..
-C     .. Local Arrays ..
-      REAL f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype)
-c
-c--->    open the file with forces
-c
-      OPEN (43,file='forces.dat',status='unknown',form='formatted')
-      REWIND 43
-
-      n=3*ntype
-      istep0=0
-
-      !The might be a line like
-      !maxhist=3
-      !at the beginning of the file
-      READ(43,'(a)',end=100) testchar
-      REWIND 43
-      IF (testchar =='m') THEN
-         READ(43,"(8x,i5)") maxhist
-         DO
-            READ(43,*,END = 99) !read the line containing the energy
-            istep0 = istep0+1
-            DO j = 1,ntype
-               READ(43,*,err= 100, end=100)
-            ENDDO
-         ENDDO
- 99      REWIND 43
-         read(43,*) !skip maxhist line
-         !skip all but last maxhist entries
-         DO i = 1, MAX(istep0-maxhist,0)
-            read(43,*,end=100)
-            DO j = 1,ntype
-               READ(43,*,err= 100, end=100)
-            ENDDO
-         ENDDO
-      ENDIF
-      istep0 = 0
-
-      DO !  read from file until end of data
-         l_skip = .FALSE.
-         READ(43,'(33x,l1)',END = 101,err = 101) l_skip
- 101     READ (43,'(3f16.9,3f14.9)',END = 100,err = 100) 
-     >        ((xnew(i+3*j),i = 1,3),(f(i+3*j),i = 1,3),j = 0,ntype-1) 
-         IF (l_skip) CYCLE !skip this entry
-
-         istep0 = istep0+1
-
-         IF (istep0.EQ.1) THEN
-!------------------------------------------------------------
-! ===>    first step
-!     
-!---  >       set h to identity matrix
-            h = 0.0
-            DO j = 1,n
-               h(j,j) = 1.0
-            ENDDO
-
-c------------------------------------------------------------
-         ELSE
-! ===>    k-th step
-            
-!---  >       determine p and shift x as needed
-            p(:) = xnew(:)-xold(:)
-!---  >       update the change in gradients
-            y(:) = y(:)-f(:)
-!---  >       get necessary inner products and H|y>
-            py = dot_product(p,y)
-            yy = 0.0
-            DO i = 1,n
-               s = 0.0
-               DO j = 1,n
-                  s = s+y(j)*h(j,i)
-               ENDDO
-               v(i) = s
-               yy = yy+y(i)*s
-            ENDDO
-!---  >       check that update will leave h positive definite,
-!---  >       if not, restart with diagonal matrix
-            IF (py<=0.0) THEN
-               WRITE (6,*) 'bfgs0: <p|y> < 0, restart istep0 =',istep0
-               h = 0.0
-               DO j = 1,n
-                  h(j,j) = 1.0
-               ENDDO
-               istep0 = 1
-            ELSE
-!---  >          for second step, use oren-spedicato scaling
-!---  >          for standard scaling, use gamma = 1.
-               IF (istep0  == 2) THEN
-                  gamma = py/yy
-               ELSE
-                  gamma = 1.0
-               ENDIF
-!---  >          update h
-               DO j = 1,n
-                  DO i = 1,n
-                     h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma
-     +                    + (1.+gamma*yy/py)*p(i)*p(j)/py
-                  ENDDO
-               ENDDO
-            ENDIF
-         ENDIF
-
-!------------------------------------------------------------
-!---  >       put y in the correct form for the next step
-         DO i = 1,n
-            xold(i) = xnew(i)
-            y(i) = f(i)
-         ENDDO
-         
-      ENDDO
-         
- 100  CLOSE(43)
-      IF (istep0==0) THEN
-         WRITE (6,1000)
-      ELSE
-         WRITE (6,1001) istep0
-      ENDIF
- 1000 FORMAT(' bfgs0: No previous atomic relaxation data found'/)
- 1001 FORMAT(' bfgs0: Approximate hessian constructed using',i4,
-     +       ' steps'/)
-
-      END SUBROUTINE bfgs0
-      END MODULE m_bfgs0

force/bfgs0.f90

+MODULE m_bfgs0
+
+CONTAINS
+
+!*******************************************************************
+!     checks whether a file exists with positions and forces to
+!     use to relax structure using the BFGS quasi-Newton method.
+!     istep0 is the number of steps included on file.
+!     h,y, and xold will be saved correctly for further use.
+!*******************************************************************
+SUBROUTINE bfgs0(ntype,istep0,xold,y,h)
+      
+   IMPLICIT NONE
+
+   ! Scalar Arguments
+   INTEGER, INTENT (IN) :: ntype
+   INTEGER, INTENT (OUT):: istep0
+
+   ! Array Arguments
+   REAL,    INTENT (OUT):: y(3*ntype),xold(3*ntype)
+   REAL,    INTENT (OUT):: h(3*ntype,3*ntype)
+
+   ! Local Scalars
+   INTEGER n,is,i,j
+   REAL yy,py,s,gamma
+   CHARACTER :: testchar
+   INTEGER   :: maxhist
+   LOGICAL   :: l_skip
+
+   ! Local Arrays
+   REAL f(3*ntype),p(3*ntype),v(3*ntype),xnew(3*ntype)
+
+   ! open the file with forces
+
+   OPEN (43,file='forces.dat',status='unknown',form='formatted')
+   REWIND 43
+
+   n=3*ntype
+   istep0=0
+
+   ! There might be a line like