relaxation.F90 5.93 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.
!--------------------------------------------------------------------------------

Daniel Wortmann's avatar
Daniel Wortmann committed
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 ...

Daniel Wortmann's avatar
Daniel Wortmann committed
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
Daniel Wortmann's avatar
Daniel Wortmann committed
19
    USE m_types
20
    USE m_relaxio
21
    USE m_mixing_history
22 23 24
#ifdef CPP_MPI
    INCLUDE 'mpif.h'
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
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

Daniel Wortmann's avatar
Daniel Wortmann committed
32 33
    REAL,ALLOCATABLE :: pos(:,:,:),force(:,:,:),energies(:)
    REAL,ALLOCATABLE :: displace(:,:),old_displace(:,:)
34 35 36
    INTEGER          :: n,ierr
    LOGICAL          :: l_conv

Daniel Wortmann's avatar
Daniel Wortmann committed
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)


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

Daniel Wortmann's avatar
Daniel Wortmann committed
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


Daniel Wortmann's avatar
Daniel Wortmann committed
100 101
    displace = alpha*force(:,:,SIZE(force,3))
  END SUBROUTINE simple_step
102

Daniel Wortmann's avatar
Daniel Wortmann committed
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(:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
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))
Daniel Wortmann's avatar
Daniel Wortmann committed
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)
Daniel Wortmann's avatar
Daniel Wortmann committed
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)
Daniel Wortmann's avatar
Daniel Wortmann committed
134
       v = MATMUL(y,h)
135
       yy = DOT_PRODUCT(y,v)
Daniel Wortmann's avatar
Daniel Wortmann committed
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))
Daniel Wortmann's avatar
Daniel Wortmann committed
162
  END SUBROUTINE simple_bfgs
163

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

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

Daniel Wortmann's avatar
Daniel Wortmann committed
173
    n_old = SIZE(pos,3)-1
174

Daniel Wortmann's avatar
Daniel Wortmann committed
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