LCOV - code coverage report
Current view: top level - force - relaxation.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 32 77 41.6 %
Date: 2019-09-08 04:53:50 Functions: 2 4 50.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_relaxation
       8             :   USE m_judft
       9             :   IMPLICIT NONE
      10             :   PRIVATE
      11             :   PUBLIC relaxation !This is the interface. Below there are internal subroutines for bfgs, simple mixing, CG ...
      12             : 
      13             : CONTAINS
      14           2 :   SUBROUTINE relaxation(mpi,input,atoms,cell,sym,force_new,energies_new)
      15             :     !This routine uses the current force,energies and atomic positions to 
      16             :     !generate a displacement in a relaxation step. 
      17             :     !The history is taken into account by read_relax from m_relaxio
      18             :     !After generating new positions the code stops
      19             :     USE m_types
      20             :     USE m_relaxio
      21             :     USE m_mixing_history
      22             : #ifdef CPP_MPI
      23             :     INCLUDE 'mpif.h'
      24             : #endif
      25             :     TYPE(t_mpi),INTENT(IN)   :: mpi
      26             :     TYPE(t_input),INTENT(IN) :: input
      27             :     TYPE(t_atoms),INTENT(IN) :: atoms
      28             :     TYPE(t_sym),INTENT(IN)   :: sym
      29             :     TYPE(t_cell),INTENT(IN)  :: cell
      30             :     REAL,INTENT(in)          :: force_new(:,:),energies_new !data for this iteration
      31             : 
      32           4 :     REAL,ALLOCATABLE :: pos(:,:,:),force(:,:,:),energies(:)
      33           4 :     REAL,ALLOCATABLE :: displace(:,:),old_displace(:,:)
      34             :     INTEGER          :: n,ierr
      35             :     LOGICAL          :: l_conv
      36             : 
      37           2 :     IF (mpi%irank==0) THEN
      38           1 :        ALLOCATE(pos(3,atoms%ntype,1)); 
      39           3 :        DO n=1,atoms%ntype
      40           3 :           pos(:,n,1)=atoms%pos(:,SUM(atoms%neq(:n-1))+1)
      41             :        END DO
      42           1 :        ALLOCATE(force(3,atoms%ntype,1)); force(:,:,1)=force_new
      43           1 :        ALLOCATE(energies(1));energies(1)=energies_new
      44           1 :        ALLOCATE(displace(3,atoms%ntype),old_displace(3,atoms%ntype))
      45             : 
      46             :        !Remove force components that are not selected for relaxation
      47           3 :        DO n=1,atoms%ntype
      48           3 :           IF (atoms%l_geo(n)) THEN
      49           2 :              force(:,n,1)=force(:,n,1)*REAL(atoms%relax(:,n))
      50             :           ELSE
      51           0 :              force(:,n,1)=0.0
      52             :           ENDIF
      53             :        ENDDO
      54             :        
      55             :        ! add history 
      56           1 :        CALL read_relax(pos,force,energies)
      57             : 
      58             :        !determine new positions
      59           1 :        IF (SIZE(energies)==1.OR.input%forcemix==0) THEN
      60             :           !no history present simple step
      61             :           ! choose a reasonable first guess for scaling
      62             :           ! this choice is based on a Debye temperature of 330K;
      63             :           ! modify as needed
      64             :           !alpha = (250.0/(MAXVAL(atoms%zatom)*input%xa))*((330./input%thetad)**2)
      65           1 :           CALL simple_step(input%forcealpha,0.25,force,displace)
      66           0 :        ELSE IF (input%forcemix==1) THEN
      67           0 :           CALL simple_cg(pos,force,displace)
      68           0 :        ELSE IF (input%forcemix==2) THEN
      69           0 :           CALL simple_bfgs(pos,force,displace)
      70             :        ELSE
      71           0 :           CALL juDFT_error('unkown mixing scheme for forces', calledby='relaxation')
      72             :        END IF
      73             : 
      74             :        !Check for convergence of forces/displacements
      75           1 :        l_conv=.TRUE.
      76           3 :        DO n=1,atoms%ntype
      77           2 :           IF (DOT_PRODUCT(force(:,n,SIZE(force,3)),force(:,n,SIZE(force,3)))>input%epsforce**2) l_conv=.FALSE.
      78           3 :           IF (DOT_PRODUCT(displace(:,n),displace(:,n))>input%epsdisp**2) l_conv=.FALSE.
      79             :        ENDDO
      80             : 
      81             :        !New displacements relative to positions in inp.xml
      82           1 :        CALL read_displacements(atoms,old_displace)
      83           1 :        displace=displace+old_displace
      84             : 
      85             :        !Write file
      86           1 :        CALL write_relax(pos,force,energies,displace)
      87             : 
      88             : 
      89             :     ENDIF
      90             : #ifdef CPP_MPI
      91           2 :     CALL MPI_BCAST(l_conv,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
      92             : #endif
      93           2 :     IF (l_conv) THEN
      94           0 :        CALL judft_end("Structural relaxation: Done",mpi%irank)
      95             :     ELSE
      96           2 :        CALL mixing_history_reset(mpi)
      97           2 :        CALL judft_end("Structural relaxation: new displacements generated",mpi%irank)
      98             :     END IF
      99           0 :   END SUBROUTINE relaxation
     100             : 
     101             : 
     102             : 
     103           1 :   SUBROUTINE simple_step(alpha,maxdisp,force,displace)
     104             :     !-----------------------------------------------
     105             :     IMPLICIT NONE
     106             :     REAL,INTENT(in)  :: alpha,maxdisp
     107             :     REAL,INTENT(in)  :: force(:,:,:)
     108             :     REAL,INTENT(OUT) :: displace(:,:)
     109             : 
     110             :     real :: corr
     111             :     
     112           1 :     displace = alpha*force(:,:,SIZE(force,3))
     113           3 :     corr=maxdisp/maxval(abs(displace))
     114           1 :     if (corr<1.0) displace = corr*alpha*force(:,:,size(force,3))
     115             :     
     116           1 :   END SUBROUTINE simple_step
     117             : 
     118           0 :   SUBROUTINE simple_bfgs(pos,force,shift)
     119             :     !-----------------------------------------------
     120             :     !  Simple BFGS method to calculate shift out of old positions and forces
     121             :     !-----------------------------------------------
     122             :     IMPLICIT NONE
     123             :     REAL,INTENT(in)  :: pos(:,:,:),force(:,:,:)
     124             :     REAL,INTENT(OUT) :: shift(:,:)
     125             : 
     126             :     INTEGER         :: n,i,j,hist_length,n_force
     127           0 :     REAL,ALLOCATABLE:: h(:,:)
     128           0 :     REAL,ALLOCATABLE:: p(:),y(:),v(:)
     129             :     REAL            :: py,yy,gamma
     130             : 
     131           0 :     n_force=3*SIZE(pos,2)
     132           0 :     ALLOCATE(h(n_force,n_force))
     133           0 :     ALLOCATE(p(n_force),y(n_force),v(n_force))
     134             : 
     135             :     !calculate approx. Hessian
     136             :     !initialize H
     137           0 :     h = 0.0
     138           0 :     DO n = 1,n_force
     139           0 :        h(n,n) = 1.0
     140             :     ENDDO
     141             :     !loop over all iterations (including current)
     142           0 :     hist_length=SIZE(pos,3)
     143           0 :     DO n = 2,hist_length
     144             :        ! differences
     145           0 :        p(:) = RESHAPE(pos(:,:,n)-pos(:,:,n-1),(/SIZE(p)/))
     146           0 :        y(:) = RESHAPE(force(:,:,n-1)-force(:,:,n),(/SIZE(p)/))
     147             :        ! get necessary inner products and H|y>
     148           0 :        py = DOT_PRODUCT(p,y)
     149           0 :        v = MATMUL(y,h)
     150           0 :        yy = DOT_PRODUCT(y,v)
     151             :        !check that update will leave h positive definite;
     152           0 :        IF (py <= 0.0) THEN
     153           0 :           WRITE (6,*) '  bfgs: <p|y> < 0'
     154           0 :           WRITE (6,*) '  check convergence of forces'
     155             :           !Starting over with initial hessian
     156           0 :           h = 0.0
     157           0 :           DO j = 1,n_force
     158           0 :              h(j,j) = 1.0
     159             :           ENDDO
     160             :           CYCLE 
     161             :        ELSE
     162             :           !update h
     163           0 :           IF (n == 2) THEN
     164           0 :              gamma = py/yy
     165             :           ELSE
     166             :              gamma = 1.0
     167             :           ENDIF
     168           0 :           DO j = 1,n_force
     169           0 :              DO i = 1,n_force
     170           0 :                 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
     171             :              ENDDO
     172             :           ENDDO
     173             :        ENDIF
     174             :     ENDDO
     175           0 :     y(:) = RESHAPE(force(:,:,hist_length),(/SIZE(p)/))
     176           0 :     shift = RESHAPE(MATMUL(y,h),SHAPE(shift))
     177           0 :   END SUBROUTINE simple_bfgs
     178             : 
     179           0 :   SUBROUTINE simple_cg(pos,force,shift)
     180             :     !-----------------------------------------------
     181             :     IMPLICIT NONE
     182             :     REAL,INTENT(in)  :: pos(:,:,:),force(:,:,:)
     183             :     REAL,INTENT(OUT) :: shift(:,:)
     184             : 
     185           0 :     REAL                :: f1(3,SIZE(pos,2)),f2(3,SIZE(pos,2))
     186             :     INTEGER             :: n_old
     187             : 
     188           0 :     n_old = SIZE(pos,3)-1
     189             : 
     190           0 :     f1 = (force(:,:,n_old+1)-force(:,:,n_old))/(pos(:,:,n_old+1)-pos(:,:,n_old))
     191           0 :     f2 = force(:,:,n_old+1)-f1*pos(:,:,n_old+1)
     192           0 :     shift = -1.*f2/f1-force(:,:,n_old+1)
     193           0 :   END SUBROUTINE simple_cg
     194             : END MODULE m_relaxation

Generated by: LCOV version 1.13