relaxation.F90 6.47 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
    IF (mpi%irank==0) THEN
       ALLOCATE(pos(3,atoms%ntype,1)); 
       DO n=1,atoms%ntype
40
          pos(:,n,1)=atoms%pos(:,SUM(atoms%neq(:n-1))+1)
41 42 43 44
       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
       !Remove force components that are not selected for relaxation
       DO n=1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
             force(:,n,1)=force(:,n,1)*REAL(atoms%relax(:,n))
          ELSE
             force(:,n,1)=0.0
          ENDIF
       ENDDO
       
55 56 57 58 59 60 61 62 63 64
       ! 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)
65
          CALL simple_step(input%forcealpha,0.25,force,displace)
66
       ELSE IF (input%forcemix==1) THEN
67
          CALL simple_cg(pos,force,displace)
68
       ELSE IF (input%forcemix==2) THEN
69
          CALL simple_bfgs(pos,force,displace)
70 71 72
       ELSE
          CALL juDFT_error('unkown mixing scheme for forces', calledby='relaxation')
       END IF
73 74 75 76 77

       !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.
78
          IF (DOT_PRODUCT(displace(:,n),displace(:,n))>input%epsdisp**2) l_conv=.FALSE.
79 80 81 82 83 84 85 86 87 88
       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)


89
    ENDIF
90
#ifdef CPP_MPI
91
    CALL MPI_BCAST(l_conv,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
92 93
#endif
    IF (l_conv) THEN
94
       CALL judft_end("Structual relaxation: Done",mpi%irank)
95
    ELSE
96
       CALL mixing_history_reset(mpi)
97
       CALL judft_end("Structual relaxation: new displacements generated",mpi%irank)
98
    END IF
99
  END SUBROUTINE relaxation
100

101 102


103
  SUBROUTINE simple_step(alpha,maxdisp,force,displace)
104 105
    !-----------------------------------------------
    IMPLICIT NONE
106
    REAL,INTENT(in)  :: alpha,maxdisp
107 108
    REAL,INTENT(in)  :: force(:,:,:)
    REAL,INTENT(OUT) :: displace(:,:)
109

110 111
    real :: corr
    
112
    displace = alpha*force(:,:,SIZE(force,3))
113 114 115
    corr=maxdisp/maxval(abs(displace))
    if (corr<1.0) displace = corr*alpha*force(:,:,size(force,3))
    
116
  END SUBROUTINE simple_step
117

118 119 120 121 122 123
  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(:,:,:)
124
    REAL,INTENT(OUT) :: shift(:,:)
125 126 127 128 129 130

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

131 132 133
    n_force=3*SIZE(pos,2)
    ALLOCATE(h(n_force,n_force))
    ALLOCATE(p(n_force),y(n_force),v(n_force))
134 135 136 137 138 139 140 141

    !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)
142
    hist_length=SIZE(pos,3)
143 144 145 146 147
    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>
148
       py = DOT_PRODUCT(p,y)
149
       v = MATMUL(y,h)
150
       yy = DOT_PRODUCT(y,v)
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
       !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)/))
176
    shift = RESHAPE(MATMUL(y,h),SHAPE(shift))
177
  END SUBROUTINE simple_bfgs
178

179 180 181
  SUBROUTINE simple_cg(pos,force,shift)
    !-----------------------------------------------
    IMPLICIT NONE
182 183
    REAL,INTENT(in)  :: pos(:,:,:),force(:,:,:)
    REAL,INTENT(OUT) :: shift(:,:)
184 185 186

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

188
    n_old = SIZE(pos,3)-1
189

190 191 192 193 194
    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