bfgs0.f90 3.8 KB
Newer Older
Gregor Michalicek's avatar
Gregor Michalicek committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
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
   ! 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)
         END DO
      END DO
      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)
         END DO
      END DO
   END IF
   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
         END DO

      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)
            END DO
            v(i) = s
            yy = yy+y(i)*s
         END DO
         ! 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
            END DO
            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
            END IF
            ! 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
               END DO
            END DO
         END IF
      END IF

      ! put y in the correct form for the next step
      DO i = 1,n
         xold(i) = xnew(i)
         y(i) = f(i)
      END DO
         
   END DO
         
   100 CLOSE(43)

   IF (istep0==0) THEN
      WRITE (6,1000)
   ELSE
      WRITE (6,1001) istep0
   END IF

   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