relaxation.F90 5.91 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10
MODULE m_relaxation
  USE m_judft
  IMPLICIT NONE
  PRIVATE
11 12
  PUBLIC relaxation !This is the interface. Below there are internal subroutines for bfgs, simple mixing, CG ...

13 14
CONTAINS
  SUBROUTINE relaxation(mpi,input,atoms,cell,sym,force_new,energies_new)
15 16 17 18
    !This routine uses the current force,energies and atomic positions to 
    !generate a displacement in a relaxation step. 
    !The history is taken into account by read_relax from m_relaxio
    !After generating new positions the code stops
19
    USE m_types
20
    USE m_relaxio
21
    USE m_mixing_history
22 23 24
#ifdef CPP_MPI
    INCLUDE 'mpif.h'
#endif
25 26 27 28 29
    TYPE(t_mpi),INTENT(IN)   :: mpi
    TYPE(t_input),INTENT(IN) :: input
    TYPE(t_atoms),INTENT(IN) :: atoms
    TYPE(t_sym),INTENT(IN)   :: sym
    TYPE(t_cell),INTENT(IN)  :: cell
30 31
    REAL,INTENT(in)          :: force_new(:,:),energies_new !data for this iteration

32 33
    REAL,ALLOCATABLE :: pos(:,:,:),force(:,:,:),energies(:)
    REAL,ALLOCATABLE :: displace(:,:),old_displace(:,:)
34 35 36
    INTEGER          :: n,ierr
    LOGICAL          :: l_conv

37 38 39 40 41 42 43 44
    IF (mpi%irank==0) THEN
       ALLOCATE(pos(3,atoms%ntype,1)); 
       DO n=1,atoms%ntype
          pos(:,n,1)=atoms%taual(:,SUM(atoms%neq(:n-1))+1)
       END DO
       ALLOCATE(force(3,atoms%ntype,1)); force(:,:,1)=force_new
       ALLOCATE(energies(1));energies(1)=energies_new
       ALLOCATE(displace(3,atoms%ntype),old_displace(3,atoms%ntype))
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

       ! add history 
       CALL read_relax(pos,force,energies)

       !determine new positions
       IF (SIZE(energies)==1.OR.input%forcemix==0) THEN
          !no history present simple step
          ! choose a reasonable first guess for scaling
          ! this choice is based on a Debye temperature of 330K;
          ! modify as needed
          !alpha = (250.0/(MAXVAL(atoms%zatom)*input%xa))*((330./input%thetad)**2)
          CALL simple_step(input%forcealpha,force,displace)
       ELSEIF (input%forcemix==1) THEN
          CALL simple_cg(pos,force,displace)
       ELSE
          CALL simple_bfgs(pos,force,displace)
       ENDIF

       !Check for convergence of forces/displacements
       l_conv=.TRUE.
       DO n=1,atoms%ntype
          IF (DOT_PRODUCT(force(:,n,SIZE(force,3)),force(:,n,SIZE(force,3)))>input%epsforce**2) l_conv=.FALSE.
          IF (DOT_PRODUCT(displace(:,n),displace(:,n))>input%epsforce**2) l_conv=.FALSE.
       ENDDO

       !New displacements relative to positions in inp.xml
       CALL read_displacements(atoms,old_displace)
       displace=displace+old_displace

       !Write file
       CALL write_relax(pos,force,energies,displace)


78
    ENDIF
79
#ifdef CPP_MPI
80
    CALL MPI_BCAST(l_conv,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
81 82 83 84
#endif
    IF (l_conv) THEN
       CALL judft_end("Structual relaxation: Done",0)
    ELSE
85
       CALL mixing_history_reset(mpi)
86 87
       CALL judft_end("Structual relaxation: new displacements generated",0)
    END IF
88
  END SUBROUTINE relaxation
89

90 91 92 93 94 95 96 97


  SUBROUTINE simple_step(alpha,force,displace)
    !-----------------------------------------------
    IMPLICIT NONE
    REAL,INTENT(in)  :: alpha
    REAL,INTENT(in)  :: force(:,:,:)
    REAL,INTENT(OUT) :: displace(:,:)
98 99


100 101
    displace = alpha*force(:,:,SIZE(force,3))
  END SUBROUTINE simple_step
102

103 104 105 106 107 108
  SUBROUTINE simple_bfgs(pos,force,shift)
    !-----------------------------------------------
    !  Simple BFGS method to calculate shift out of old positions and forces
    !-----------------------------------------------
    IMPLICIT NONE
    REAL,INTENT(in)  :: pos(:,:,:),force(:,:,:)
109
    REAL,INTENT(OUT) :: shift(:,:)
110 111 112 113 114 115

    INTEGER         :: n,i,j,hist_length,n_force
    REAL,ALLOCATABLE:: h(:,:)
    REAL,ALLOCATABLE:: p(:),y(:),v(:)
    REAL            :: py,yy,gamma

116 117 118
    n_force=3*SIZE(pos,2)
    ALLOCATE(h(n_force,n_force))
    ALLOCATE(p(n_force),y(n_force),v(n_force))
119 120 121 122 123 124 125 126

    !calculate approx. Hessian
    !initialize H
    h = 0.0
    DO n = 1,n_force
       h(n,n) = 1.0
    ENDDO
    !loop over all iterations (including current)
127
    hist_length=SIZE(pos,3)
128 129 130 131 132
    DO n = 2,hist_length
       ! differences
       p(:) = RESHAPE(pos(:,:,n)-pos(:,:,n-1),(/SIZE(p)/))
       y(:) = RESHAPE(force(:,:,n)-force(:,:,n-1),(/SIZE(p)/))
       ! get necessary inner products and H|y>
133
       py = DOT_PRODUCT(p,y)
134
       v = MATMUL(y,h)
135
       yy = DOT_PRODUCT(y,v)
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
       !check that update will leave h positive definite;
       IF (py <= 0.0) THEN
          WRITE (6,*) '  bfgs: <p|y> < 0'
          WRITE (6,*) '  check convergence of forces'
          !Starting over with initial hessian
          h = 0.0
          DO j = 1,n_force
             h(j,j) = 1.0
          ENDDO
          CYCLE 
       ELSE
          !update h
          IF (n == 2) THEN
             gamma = py/yy
          ELSE
             gamma = 1.0
          ENDIF
          DO j = 1,n_force
             DO i = 1,n_force
                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
    ENDDO
    y(:) = RESHAPE(force(:,:,hist_length),(/SIZE(p)/))
161
    shift = RESHAPE(MATMUL(y,h),SHAPE(shift))
162
  END SUBROUTINE simple_bfgs
163

164 165 166
  SUBROUTINE simple_cg(pos,force,shift)
    !-----------------------------------------------
    IMPLICIT NONE
167 168
    REAL,INTENT(in)  :: pos(:,:,:),force(:,:,:)
    REAL,INTENT(OUT) :: shift(:,:)
169 170 171

    REAL                :: f1(3,SIZE(pos,2)),f2(3,SIZE(pos,2))
    INTEGER             :: n_old
172

173
    n_old = SIZE(pos,3)-1
174

175 176 177 178 179
    f1 = (force(:,:,n_old+1)-force(:,:,n_old))/(pos(:,:,n_old+1)-pos(:,:,n_old))
    f2 = force(:,:,n_old+1)-f1*pos(:,:,n_old+1)
    shift = -1.*f2/f1-force(:,:,n_old+1)
  END SUBROUTINE simple_cg
END MODULE m_relaxation