bfgs.f90 5.08 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 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
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