relaxation.F90 7.01 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
    !This routine uses the current force,energies and atomic positions to
    !generate a displacement in a relaxation step.
17 18
    !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
    USE m_types_xml
23 24 25
#ifdef CPP_MPI
    INCLUDE 'mpif.h'
#endif
26 27 28 29 30
    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
31 32
    REAL,INTENT(in)          :: force_new(:,:),energies_new !data for this iteration

33 34
    REAL,ALLOCATABLE :: pos(:,:,:),force(:,:,:),energies(:)
    REAL,ALLOCATABLE :: displace(:,:),old_displace(:,:)
35 36
    INTEGER          :: n,ierr
    LOGICAL          :: l_conv
37 38 39
    !to calculate the current displacement
    Type(t_xml)  :: xml
    TYPE(t_atoms):: atoms_non_displaced
40

41
    IF (mpi%irank==0) THEN
42 43
       call xml%init()
       ALLOCATE(pos(3,atoms%ntype,1));
44
       DO n=1,atoms%ntype
45
          pos(:,n,1)=atoms%pos(:,SUM(atoms%neq(:n-1))+1)
46 47 48 49
       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))
50

51 52 53 54 55 56 57 58
       !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
59 60

       ! add history
61 62 63 64 65 66 67 68 69
       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)
70
          CALL simple_step(input%forcealpha,0.25,force,displace)
71
       ELSE IF (input%forcemix==1) THEN
72
          CALL simple_cg(pos,force,displace)
73
       ELSE IF (input%forcemix==2) THEN
74
          CALL simple_bfgs(pos,force,displace)
75
       ELSE
76
          CALL juDFT_error('Unknown mixing scheme for forces', calledby='relaxation')
77
       END IF
78 79 80 81 82

       !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.
83
          IF (DOT_PRODUCT(displace(:,n),displace(:,n))>input%epsdisp**2) l_conv=.FALSE.
84 85 86
       ENDDO

       !New displacements relative to positions in inp.xml
87 88 89 90 91
       !CALL read_displacements(atoms,old_displace)
       call atoms_non_displaced%read_xml(xml)
       call xml%freeResources()
       old_displace=atoms%taual-atoms_non_displaced%taual

92 93 94 95 96 97
       displace=displace+old_displace

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


98
    ENDIF
99
#ifdef CPP_MPI
100
    CALL MPI_BCAST(l_conv,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
101 102
#endif
    IF (l_conv) THEN
103
       CALL judft_end("Structural relaxation: Done",mpi%irank)
104
    ELSE
105
       CALL mixing_history_reset(mpi)
106
       CALL judft_end("Structural relaxation: new displacements generated",mpi%irank)
107
    END IF
108
  END SUBROUTINE relaxation
109

110 111


112
  SUBROUTINE simple_step(alpha,maxdisp,force,displace)
113 114
    !-----------------------------------------------
    IMPLICIT NONE
115
    REAL,INTENT(in)  :: alpha,maxdisp
116 117
    REAL,INTENT(in)  :: force(:,:,:)
    REAL,INTENT(OUT) :: displace(:,:)
118

119
    real :: corr
120

121
    displace = alpha*force(:,:,SIZE(force,3))
122 123
    corr=maxdisp/maxval(abs(displace))
    if (corr<1.0) displace = corr*alpha*force(:,:,size(force,3))
124

125
  END SUBROUTINE simple_step
126

127 128 129 130 131 132
  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(:,:,:)
133
    REAL,INTENT(OUT) :: shift(:,:)
134 135 136 137 138 139

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

140 141 142
    n_force=3*SIZE(pos,2)
    ALLOCATE(h(n_force,n_force))
    ALLOCATE(p(n_force),y(n_force),v(n_force))
143 144 145 146 147 148 149 150

    !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)
151
    hist_length=SIZE(pos,3)
152 153 154
    DO n = 2,hist_length
       ! differences
       p(:) = RESHAPE(pos(:,:,n)-pos(:,:,n-1),(/SIZE(p)/))
155
       y(:) = RESHAPE(force(:,:,n-1)-force(:,:,n),(/SIZE(p)/))
156
       ! get necessary inner products and H|y>
157
       py = DOT_PRODUCT(p,y)
158
       v = MATMUL(y,h)
159
       yy = DOT_PRODUCT(y,v)
160 161 162 163 164 165 166 167 168
       !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
169
          CYCLE
170 171 172 173 174 175 176 177 178 179 180 181 182 183 184
       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)/))
185
    shift = RESHAPE(MATMUL(y,h),SHAPE(shift))
186
  END SUBROUTINE simple_bfgs
187

188 189 190
  SUBROUTINE simple_cg(pos,force,shift)
    !-----------------------------------------------
    IMPLICIT NONE
191 192
    REAL,INTENT(in)  :: pos(:,:,:),force(:,:,:)
    REAL,INTENT(OUT) :: shift(:,:)
193 194

    REAL                :: f1(3,SIZE(pos,2)),f2(3,SIZE(pos,2))
195 196 197 198 199
    REAL                :: dist(3,SIZE(pos,2))
    REAL                :: eps
    INTEGER             :: n_old, i, j

    eps = 1.0e-9
200

201
    n_old = SIZE(pos,3)-1
202

203 204 205
    dist(:,:) = pos(:,:,n_old+1)-pos(:,:,n_old)
    DO i = 1, SIZE(pos,2)
       DO j = 1, 3
206
          IF(ABS(dist(j,i)).LT.eps) dist(j,i) = eps ! To avoid calculation of 0.0/0.0 below.
207 208 209 210
       END DO
    END DO

    f1 = (force(:,:,n_old+1)-force(:,:,n_old))/dist
211 212 213 214
    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