force_w.F90 4.42 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
!--------------------------------------------------------------------------------
! Copyright (c) 2019 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.
!--------------------------------------------------------------------------------
MODULE m_forcew
  ! ************************************************************
  ! Printing force components
  ! ************************************************************
CONTAINS
  SUBROUTINE force_w(mpi,input,atoms,sym,results,cell,oneD,vacuum)
    USE m_types
    USE m_xmlOutput
    USE m_relaxation
    IMPLICIT NONE
    TYPE(t_mpi),INTENT(IN)       :: mpi
    TYPE(t_results),INTENT(INOUT):: results
    TYPE(t_oneD),INTENT(IN)      :: oneD
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_sym),INTENT(IN)       :: sym
    TYPE(t_cell),INTENT(IN)      :: cell
    TYPE(t_atoms),INTENT(IN)     :: atoms
    TYPE(t_vacuum),INTENT(IN)    :: vacuum
    !     ..
    !     .. Local Scalars ..
    REAL,PARAMETER:: zero=0.0
    REAL sum
    INTEGER i,jsp,n,nat1,ierr
    REAL eps_force
    LOGICAL :: l_new,l_relax
    !     ..
    !     .. Local Arrays ..
    REAL forcetot(3,atoms%ntype)
    CHARACTER(LEN=20) :: attributes(7)
35 36 37
#ifdef CPP_MPI
    include 'mpif.h'
#endif
38 39 40 41 42 43 44 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
    !
    !     write spin-dependent forces
    !
    IF (mpi%irank==0) THEN
       nat1 = 1
       DO n = 1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
             IF (input%jspins.EQ.2) THEN
                DO jsp = 1,input%jspins
                   WRITE (6,FMT=8000) jsp,n, (atoms%pos(i,nat1),i=1,3),&
                        &           (results%force(i,n,jsp),i=1,3)
                END DO
             END IF
8000         FORMAT ('SPIN-',i1,1x,'FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,&
                  &          3x,'Y=',f7.3,3x,'Z=',f7.3,5x,' FX_SP =',f9.6,' FY_SP =',&
                  &          f9.6,' FZ_SP =',f9.6)
          ENDIF
          nat1 = nat1 + atoms%neq(n)
       END DO
       !
       !     write total forces
       !
       WRITE  (6,8005)
8005   FORMAT (/,' ***** TOTAL FORCES ON ATOMS ***** ',/)
       IF (input%l_f) CALL openXMLElement('totalForcesOnRepresentativeAtoms',(/'units'/),(/'Htr/bohr'/))
       nat1 = 1
       DO n = 1,atoms%ntype
          IF (atoms%l_geo(n)) THEN
             DO i = 1,3
                forcetot(i,n) = zero
             END DO
             DO jsp = 1,input%jspins
                DO i = 1,3
                   forcetot(i,n) = forcetot(i,n) + results%force(i,n,jsp)
                END DO
             END DO
74

75 76 77 78 79
             WRITE (6,FMT=8010) n, (atoms%pos(i,nat1),i=1,3),&
                  &        (forcetot(i,n),i=1,3)
8010         FORMAT (' TOTAL FORCE FOR ATOM TYPE=',i3,2x,'X=',f7.3,3x,'Y=',&
                  &              f7.3,3x,'Z=',f7.3,/,22x,' FX_TOT=',f9.6,&
                  &              ' FY_TOT=',f9.6,' FZ_TOT=',f9.6)
80

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
             WRITE(attributes(1),'(i0)') n
             WRITE(attributes(2),'(f12.6)') atoms%pos(1,nat1)
             WRITE(attributes(3),'(f12.6)') atoms%pos(2,nat1)
             WRITE(attributes(4),'(f12.6)') atoms%pos(3,nat1)
             WRITE(attributes(5),'(f12.8)') forcetot(1,n)
             WRITE(attributes(6),'(f12.8)') forcetot(2,n)
             WRITE(attributes(7),'(f12.8)') forcetot(3,n)
             IF (input%l_f) THEN
                CALL writeXMLElementFormPoly('forceTotal',(/'atomType','x       ','y       ','z       ',&
                     'F_x     ','F_y     ','F_z     '/),attributes,&
                     RESHAPE((/8,1,1,1,3,3,3,6,12,12,12,12,12,12/),(/7,2/)))
             END IF
          END IF
          nat1 = nat1 + atoms%neq(n)
       END DO
       IF (input%l_f) CALL closeXMLElement('totalForcesOnRepresentativeAtoms')
97

98 99 100 101 102 103 104 105 106 107 108

       !Check convergence of force by comparing force with old_force
       sum=MAXVAL(ABS(forcetot - results%force_old))
       results%force_old(:,:)=forcetot !Store for next iteration
       results%force=0.0
       l_relax=sum<input%force_converged
       IF (.NOT.l_relax) THEN
          WRITE (6,8020) input%force_converged,sum
8020      FORMAT ('No new postions, force convergence required=',f8.5,'max force distance=',f8.5)
       END IF
    ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
109
#ifdef CPP_MPI
110
    CALL MPI_BCAST(l_relax,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
Daniel Wortmann's avatar
Daniel Wortmann committed
111
#endif
112 113 114
    IF (l_relax.AND.input%l_f) CALL relaxation(mpi,input,atoms,cell,sym,forcetot,results%tote)

  END SUBROUTINE force_w
Daniel Wortmann's avatar
Daniel Wortmann committed
115
END MODULE m_forcew