artn

subroutine m_artn::artn(nat, etot_eng, force, ityp, tau, order, at, if_pos, disp_code, displ_vec, lconv)

Main ARTn plugin subroutine.

Author

Matic Poberznik, Miha Gunde, Nicolas Salles, Antoine Jay

    use h_artn_precision, only: dp

    use d_artn_params, only : eigenvec, push, nevalf_max, &
                              lrestart, lrelax, llanczos, linit, lend, leigen, lperp, &
                              lmove_nextmin, lpush_over, lbackward, lpush_final, &
                              iperp, irelax, istep, iover, ifound, &
                              in_lanczos_at_min, lanczos_at_min, initpfname, etot_diff_limit,   &
                              struc_format_out, verbose, void, relx, lanc, &
                              eigenfname, fpush_factor, filout, &
                              flag_false, fill_param_step, current_disp_code

    use d_artn_data, only: save_step_data
    use d_artn_data, only: etot_step, etot_sad, etot_final, etot_init
    use d_artn_data, only: typ_init, typ_sad, typ_step
    use d_artn_data, only: tau_init, tau_sad, tau_step
    use d_artn_data, only: eigen_sad, force_step
    use d_artn_data, only: de_back, de_fwd

    use m_artn_option, only: move_nextmin, read_restart

    use m_artn_tools, only: split_field

    use m_setup_artn, only: start_guess_push, isetup

    use m_artn_report, only: write_end_report, write_comment
    use m_artn_report, only: write_struc2file
    use m_artn_report, only: write_header_report
    use m_artn_report, only: write_report, write_inter_report
    use m_artn_report, only: prev_push

    use m_block_lanczos, only: block_lanczos, ilanc, lowest_eigval
    use m_artn_tools, only: ddot, dnrm2
    ! use m_artn_tools, only: compute_delr_vec, sum_force
    !
    IMPLICIT NONE

    ! -- ARGUMENTS
    INTEGER, value,   INTENT(IN)    :: nat              !  number of atoms
    REAL(DP),         INTENT(IN)    :: etot_eng         !  total energy in current step
    REAL(DP),         INTENT(IN)    :: force(3,nat)     !  force calculated by the engine
    INTEGER,          INTENT(INOUT) :: ityp(nat)        !  atom types
    REAL(DP),         INTENT(INOUT) :: tau(3,nat)       !  atomic positions (modif after forward relax)
    INTEGER,          INTENT(IN)    :: order(nat)       !  Engine order of atom
    REAL(DP),         INTENT(IN)    :: at(3,3)          !  lattice parameters in alat units in columns:
                                                        !     at(:,1)=a, at(:,2)=b, at(:,3)=c
    INTEGER,          INTENT(IN)    :: if_pos(3,nat)    !  coordinates fixed by engine
    INTEGER,          INTENT(OUT)   :: disp_code        !  encoder of stage for move_mode
    REAL(DP),         INTENT(OUT)   :: displ_vec(3,nat) !  displacement vector communicated to move mode
    LOGICAL,          INTENT(OUT)   :: lconv            !  flag for controlling convergence

    ! -- LOCAL VARIABLES
    REAL(DP)                        :: fpara(3,nat)     ! force parallel to push/eigenvec
    REAL(DP)                        :: fperp(3,nat)     ! force parallel to push/eigenvec
    LOGICAL                         :: lforc_conv       ! flag true when forces are converged
    LOGICAL                         :: lsaddle_conv     ! flag true when saddle is reached
    integer                         :: ierr

    ! write(*,*) "enter artn with istep",istep
    !
    !*> @par The ARTn algorithm proceeds as follows:
    !*  ============================================
    !*> (1) push atoms in the direction specified by user & relax in the perpendicular direction \n
    !*> (2) use the lanczos algorithm calculate the lowest eigenvalue/eigenvec \n
    !*> (3) a negative eigenvalue, update push direction otherwise push again \n
    !*> (4) follow the lanczos direction twoard the saddle point \n
    !*> (5) push twoards adjacent minimum & initial minimum \n
    !
    ! ... Flags that controls convergence
    lconv        = .false.
    lforc_conv   = .false.
    lsaddle_conv = .false.
    !
    disp_code = void

    !! artn is already finished but called more times.
    IF( lend ) THEN
       ! write(*,*) "ARTn has already finished, RETURN"
       if( verbose > 1 ) call write_comment( trim(filout), "Enter in ARTn but already finished, RETURN")

       !! call finalize, even if not done anything, since we always need to fill the variables:
       !! disp_code, displ_vec, and lconv
       call block_finalize( .true., 0, disp_code, displ_vec )
       lconv = .true.
       return
    END IF
    !


    !! check if setup has been done or not
    if( isetup == 0 ) then
       ierr = artn_error
       call err_set(ierr, __file__,__line__,&
            msg="setup_artn has not been done before calling artn()!" )
       lconv = .true.
       goto 666
    end if



    !
    ! ...Fill variables of d_artn_params (arrays are ordered !!): needs to know engine_units
    !    The variables which are known from engine are filled:
    !        natoms, lat, etot_step, types, force_step, tau_step
    !
    ierr = fill_param_step( nat, at, order, ityp, tau, etot_eng, force )
    !! Something went wrong in filling the arrays!
    IF ( ierr /= 0 ) THEN
       call err_caller(__file__,__line__)
       lconv = .true.
       goto 666
    ENDIF


    !
    ! ... Initialize artn
    istep0: IF( istep == 0 )THEN !! ------------------------------------------------------- ISTEP = 0
       !
       lend = .false.


       !
       ! ...Start to write the output
       CALL write_header_report( )


       !!
       !! create start guess if needed
       !! NOTE: could be moved into setup?
       !!
       ierr = start_guess_push( nat, push )
       if( ierr /= 0 ) then
          call err_caller(__file__,__line__)
          lconv = .true.
          goto 666
       end if

       !
       ! ... save initial data
       ierr = save_step_data( "init" )
       if( ierr /= 0 ) then
          call err_caller(__file__,__line__)
          lconv = .true.
          goto 666
       end if


       !!
       !!==========================================
       !! restart will overwrite all params
       IF( lrestart ) THEN
          !! overwrite previously initialised things with read from restart
          !
          ! ...Signal that it is a restart
          if( verbose > 1 )call write_comment( trim(filout), "Restarted previous ARTn calculation" )
          !
          ! ...Read the FLAGS, FORCES, POSITIONS, ENERGY, ...
          ierr = read_restart( )
          IF( ierr/=0 )THEN
             call err_caller(__file__,__line__)
             lconv = .true.
             goto 666
          ENDIF
          !
          ! ...Overwirte the engine Arrays with data from restart
          tau(:,:) = tau_step(:,order(:))
          ityp(:) = typ_step(order(:))
          !
       END IF
       !!==========================================


       !
       ! ...Write the initial structure to `initpfname`
       IF (verbose>1) THEN
          call write_struc2file( "initp" )
       ENDIF

       ! open(newunit=u0,file="sscheck.xyz",status="unknown",position="append")
       ! close(u0, status="delete")

    ENDIF istep0


    !!
    !! fill delr_step (experimental)
    !!
    ! ! check alloc status of delr_vec
    ! call allocate_var( 3, nat, delr_vec, 0.0_DP )
    ! !! compute delr of current tau with tau_init
    ! call compute_delr_vec( nat, tau_step, tau_init, lat, delr_vec )
    ! call sum_force( delr_vec, nat, delr_step )
    ! open(newunit=u0,file="sscheck.xyz",status="unknown",position="append")
    ! write(u0,*) nat
    ! write(u0,*) 'Lattice="',lat,'" properties=species:I:1:pos:R:3:id:I:1:force:R:3'
    ! do i = 1, nat
    !    write(u0,*) typ_step(i), tau_step(:,i), i, delr_vec(:,i)
    ! end do
    ! close(u0)

    ! ...Split the force field in para/perp field following the push field
    CALL split_field( nat, force_step, if_pos, push, fperp, fpara )

    ! ...Write Output
    CALL write_report( etot_step, force_step, fperp, fpara, lowest_eigval, if_pos, istep, nat )

    if( istep .ne. 0 ) then
       CALL check_force_convergence( nat, force_step, if_pos, fperp, fpara, lforc_conv, lsaddle_conv )
    end if


    !
    !  the basic ARTn blocks: init, perp, eigen
    !
    IF ( linit ) THEN
       !
       ! initial displacement:
       ! if `ninit = 0`, pass directly to lanczos,
       ! else set displ_vec = push, and set `lperp=.true.`
       ierr = block_pushinit( disp_code, displ_vec )
       if( ierr /= 0 ) then
          call err_caller(__file__,__line__)
          lconv = .true.
          goto 666
       end if
       !
    ELSE IF ( lperp ) THEN
       !
       ! perpendicular relaxation:
       ! set displ_vec = fperp
       ! the lperp block flag is turned off by check_force_convergence()
       ierr = block_perprelax( nat, fperp, disp_code, displ_vec )
       !
       if( ierr /= 0 ) then
          call err_caller(__file__,__line__)
          lconv = .true.
          goto 666
       end if
       !
    ELSE IF ( leigen  )THEN
       !
       ! displacement with eigenvector
       ! set displ_vec = eigenvec*current_step_size, and set `lperp=.true.`
       ierr = block_pusheigen( disp_code, displ_vec )
       if( ierr /= 0 ) then
          call err_caller(__file__,__line__)
          lconv = .true.
          goto 666
       end if

       !
       ! Write the latest eigenvec to a file (eigenvec instead of force in arguments)
       !
       IF (verbose>1) THEN
          call write_struc2file( "latest_eigenvec" )
       ENDIF
       !
    END IF


    !
    ! The saddle point is reached -> confirmed by check_force_convergence()
    !
    IF( lsaddle_conv )THEN
       !
       ! ... write the structure to file 'outfile' = prefix_sad + nsaddle
       call write_struc2file( "saddle" )
       !
       ! save the saddle point data
       !
       ierr = save_step_data( "sad" )
       if( ierr /= 0 ) then
          call err_caller(__file__,__line__)
          lconv = .true.
          goto 666
       end if
       !
       ! switch on lpush_over block
       lpush_over = .true.
       ifound = ifound + 1
       !
       ! ...write the report
       CALL write_end_report( lpush_over, lpush_final, etot_step-etot_init )
       !
       !! If the saddle point is lower in energy
       !!  than the initial point: Mode refine
       IF ( etot_sad < etot_init .and. verbose > 1 ) THEN
          call write_comment( filout, "NOTE::E_Saddle < E_init" )
       ENDIF
       !
       ! set relevant counters to zero
    ENDIF
    !
    ! ...If saddle point is reached
    ! This block performs only the PUSH to adjacent minima after the saddle point is found
    !
    IF ( lpush_over ) THEN

       !
       !
       ! do we do a final push ?
       !
       IF ( lpush_final ) THEN
          !
          ! perform step_over
          ierr = block_pushover( disp_code, displ_vec )
          if( ierr /= 0 ) then
             call err_caller(__file__,__line__)
             lconv = .true.
             goto 666
          end if

          !
       ELSE  ! --- NO FINAL_PUSH
          !
          !! At this point the saddle point is already written by write_struct in lsaddle_conv block
          !! Here we finish the ARTn search.
          !! Preparation of the possible new ARTn search following this step.
          !! - Cleaning the flag/parameter
          !! - write in output saying no more research
          !! - return a configuration in which a new ARTn search can start
          !
          !! should be block finalzie====
          IF( verbose > 1 ) THEN
             call write_comment( filout, "NO FINAL_PUSH :: Return to the start configuration" )
          END IF

          ! ...Return to the initial comfiguration
          ! tau(:,:) = tau_init(:,order(:))
          ! ityp(:) = typ_init(order(:))

          ! put block flags to false
          call flag_false()

          ! ...Tell to the engine it is finished
          lconv = .true.

          ! ...Set the force to zero
          displ_vec(:,:) = 0.0_dp

       ENDIF
    ENDIF
    !
    ! perform a FIRE relaxation (only for reaching adjacent minima)
    !
    relax: IF ( lrelax ) THEN
       !
       ! reset force
       !
       disp_code = relx
       displ_vec = force_step
       irelax    = irelax + 1
       ilanc     = 0
       iperp     = 0
       prev_push = disp_code !! Save the previous displacement
       !
       ! The convergence is reached:
       !  - Switch the push_over or
       !  - Finish the ARTn search
       !
       IF ( lforc_conv .AND. .NOT. llanczos ) THEN
          !
          IF ( lanczos_at_min .AND. .NOT. in_lanczos_at_min ) THEN
             !
             ! ... Minimun has been found (lforc_conv =.true.)
             ! We now do lanczos to check the lowest eigenvalue
             in_lanczos_at_min = .true.
             lpush_over        = .false.
             lrelax            = .false.
             llanczos          = .true.
             disp_code         = lanc
             IF( verbose > 1 ) THEN
                call write_comment( filout, "We do a Lanczos loop at the minimum to check if lowest eivenvalue is <0" )
             END IF
             !
          ELSE
             !
             IF ( fpush_factor == 1 ) THEN
                !
                ! ... found the forward minimum!
                !   Write it to file 'outfile' = prefix_min + nmin, and return to the saddle point
                call write_struc2file( "min1" )
                !
                ! next step is relax in other direction
                disp_code = relx
                !
                ! save the min1 data
                !
                ierr = save_step_data( "min1" )
                if( ierr /= 0 ) then
                   call err_caller(__file__,__line__)
                   lconv = .true.
                   goto 666
                end if


                !
                ! ...restart from saddle point
                tau(:,:)      = tau_sad(:,order(:))
                ityp(:)       = typ_sad(order(:))
                eigenvec(:,:) = eigen_sad(:,:)
                lbackward     = .true.
                !
                ! ...Return to Push_Over Step in opposit direction
                !lsaddle = .true.
                lpush_over = .true.
                lrelax     = .false.
                !
                etot_final = etot_step
                !
                ! energy difference is saddle - current
                de_back = etot_sad - etot_final
                !
                call write_inter_report( fpush_factor, [de_back] )
                !
                ! ...reverse direction for the push_over
                fpush_factor = -1
                irelax = 0
                iover = 0
                in_lanczos_at_min = .false.
                !
             ELSE  
                !
                ! ... found the backward minimum!
                !     Write it to file 'outfile' = prefix_min + nmin
                call write_struc2file( "min2" )
                !
                ! save the min2 data
                !
                ierr = save_step_data( "min2" )
                if( ierr /= 0 ) then
                   call err_caller(__file__,__line__)
                   lconv = .true.
                   goto 666
                end if

                !
                ! ...Save the Energy difference as saddle - current
                de_fwd = etot_sad - etot_step
                !
                call write_inter_report( fpush_factor, &
                     [de_back, de_fwd, etot_init, etot_final, etot_step] )
                !
                ! signal convergence flag
                lconv = .true.
                !
             END IF
             !
          END IF
          !
       END IF
       !
    END IF relax




    !! WHAT FOR THIS BLOCK???
    !!  This should be in check_force()
    IF( etot_step - etot_init > etot_diff_limit ) then
       ierr = artn_failure
       call err_set(ierr,__file__,__line__, msg="ENERGY EXCEEDS THE LIMIT")
       lconv = .true.
       goto 666
    ENDIF

    IF( istep + 1 > nevalf_max ) then ! istep + 1 because it start at 0
       ierr = artn_failure
       call err_set(ierr, __file__,__line__, msg="NUMBER OF STEPS EXCEEDS THE LIMIT")
       lconv = .true.
       goto 666
    ENDIF




    !
    ! check if we should perform the lanczos algorithm
    !
    ! Lanczos at the end. Reason: check convergence at saddle before going into
    ! un-needed lanczos near saddle.
    !
    lanczos_: IF ( llanczos ) THEN
       !
       ierr = block_lanczos( disp_code, displ_vec, if_pos )
       !
       if( ierr /= 0 ) then
          call err_caller(__file__,__line__)
          lconv = .true.
          goto 666
       end if
       !
    ENDIF lanczos_




666 continue
    !
    !! --- Finalization Block and error check
    !
    IF( lconv )THEN
       !
       call block_finalize( lconv, ierr, disp_code, displ_vec )

       ! overwrite engine arrays
       IF( lmove_nextmin.and.ierr==0 ) then
          ! lmove_nextmin possible only if the minimization converged
          !
          ! ...Here we should load the next minimum if the user ask
          CALL move_nextmin( nat, ityp, tau, order )
       else
          !
          ! reload initial positions
          ! write(*,*) "reloading initial positions",istep
          tau(:,:) = tau_init(:,order(:))
          ityp(:) = typ_init(order(:))

          ! disp_code = RSET
          ! displ_vec = tau_init(:,:) - tau(:,:)
          ! displ_vec = displ_vec! /0.529177

          ! block
          !   integer :: i
          !   real(DP), dimension(3) :: rdum
          !   do i = 1, nat
          !      rdum = displ_vec(:,i)
          !      write(*,'(i4,1x,3f9.4,2x,f9.4)') i, rdum, norm2(rdum)
          !   end do
          ! end block

       end IF
       !
    ENDIF
    !
    ! ...Increment the ARTn-step
    istep = istep + 1

    ! save the disp_code for `extract`
    ! disp_code is passed as intent(out) argument in this routine, maybe change
    current_disp_code = disp_code
    !
Purpose

Modifies the input force to perform the ARTn algorithm

Note

d_artn_params for variables and counters that need to be stored in each step DEFINED IN: d_artn_params_mod.f90

Parameters:
  • nat[in] number of atoms

  • etot_eng[inout] total energy of the engine

  • force[in] force calculated by the engine

  • ityp[in] list of type of atoms

  • tau[inout] atomic position

  • order[in] order of atomic index in the list: force, tau, ityp

  • at[in] lattice vectors (in columns)

  • if_pos[in] list of fixed atomic dof (0 or 1)

  • disp_code[out] encoder of stage for move_mode

  • displ_vec[out] displacement vector communicated to move_mode

  • lconv[out] flag for controlling convergence

  • nat[in] [art]