Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2019 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             : MODULE m_forcew
       7             :    !-----------------------------------------------------------------------------
       8             :    ! Printing force components
       9             :    !-----------------------------------------------------------------------------
      10             : 
      11             : #ifdef CPP_MPI
      12             :    USE mpi 
      13             : #endif
      14             : 
      15             : CONTAINS
      16         668 :    SUBROUTINE force_w(fmpi,input,atoms,sym,results,cell ,vacuum)
      17             :       USE m_types
      18             :       USE m_constants
      19             :       USE m_xmlOutput
      20             :       USE m_relaxation
      21             :       USE m_rotate_forces
      22             :       IMPLICIT NONE
      23             : 
      24             :       TYPE(t_mpi),     INTENT(IN)    :: fmpi
      25             :       TYPE(t_results), INTENT(INOUT) :: results
      26             :        
      27             :       TYPE(t_input),   INTENT(IN)    :: input
      28             :       TYPE(t_sym),     INTENT(IN)    :: sym
      29             :       TYPE(t_cell),    INTENT(IN)    :: cell
      30             :       TYPE(t_atoms),   INTENT(IN)    :: atoms
      31             :       TYPE(t_vacuum),  INTENT(IN)    :: vacuum
      32             : 
      33             :       REAL maxAbsForceDist
      34             :       INTEGER i, jsp, n, nat1, ierr
      35             :       REAL eps_force
      36             :       LOGICAL :: l_new, l_forceConverged
      37         668 :       REAL forcetot(3,atoms%ntype)
      38             :       CHARACTER(LEN=20) :: attributes(7)
      39             : 
      40             :       ! Write spin-dependent forces
      41             : 
      42         668 :       IF (fmpi%irank==0) THEN
      43         919 :          DO n = 1, atoms%ntype
      44         585 :             nat1 = atoms%firstAtom(n)
      45         919 :             IF (atoms%l_geo(n)) THEN
      46         583 :                IF (input%jspins.EQ.2) THEN
      47         963 :                   DO jsp = 1,input%jspins
      48        3210 :                      WRITE (oUnit,FMT=8000) jsp, n, (atoms%pos(i,nat1),i=1,3), &
      49        6099 :                                                     (results%force(i,n,jsp),i=1,3)
      50             :                   END DO
      51             :                END IF
      52             : 
      53             : 8000           FORMAT ('SPIN-',i1,1x,'FORCE FOR ATOM TYPE=',i3,2x,&
      54             :                        'X=',f7.3,3x,'Y=',f7.3,3x,'Z=',f7.3,5x,&
      55             :                        ' FX_SP =',f9.6,' FY_SP =',f9.6,' FZ_SP =',f9.6)
      56             :             END IF
      57             :          END DO
      58             : 
      59             :          ! Write total forces
      60         334 :          WRITE  (oUnit,8005)
      61             : 
      62             : 8005     FORMAT (/,' ***** TOTAL FORCES ON ATOMS ***** ',/)
      63         392 :          IF (input%l_f) CALL openXMLElement('totalForcesOnRepresentativeAtoms',(/'units'/),(/'Htr/bohr'/))
      64        2674 :          forcetot = 0.0
      65         334 :          if (allocated(results%force_vdw)) THEN
      66           0 :             forcetot=forcetot+results%force_vdw
      67           0 :             write(oUnit,*) "vdW forces included in total force"
      68             :          endif
      69             : 
      70         919 :          DO n = 1,atoms%ntype
      71         585 :             nat1 = atoms%firstAtom(n)
      72         919 :             IF (atoms%l_geo(n)) THEN
      73        1487 :                DO jsp = 1,input%jspins
      74        4199 :                   forcetot(:,n) = forcetot(:,n) + results%force(:,n,jsp)
      75             :                END DO
      76             : 
      77             :                
      78        2915 :                WRITE (oUnit,FMT=8010) n, (atoms%pos(i,nat1),i=1,3), &
      79        3498 :                                          (forcetot(i,n),i=1,3)
      80             : 
      81             : 8010           FORMAT (' TOTAL FORCE FOR ATOM TYPE=',i3,2x,&
      82             :                        'X=',f7.3,3x,'Y=',f7.3,3x,'Z=',f7.3,/,22x,&
      83             :                        ' FX_TOT=',f9.6,' FY_TOT=',f9.6,' FZ_TOT=',f9.6)
      84             : 
      85         583 :                WRITE(attributes(1),'(i0)') n
      86         583 :                WRITE(attributes(2),'(f12.6)') atoms%pos(1,nat1)
      87         583 :                WRITE(attributes(3),'(f12.6)') atoms%pos(2,nat1)
      88         583 :                WRITE(attributes(4),'(f12.6)') atoms%pos(3,nat1)
      89         583 :                WRITE(attributes(5),'(f12.8)') forcetot(1,n)
      90         583 :                WRITE(attributes(6),'(f12.8)') forcetot(2,n)
      91         583 :                WRITE(attributes(7),'(f12.8)') forcetot(3,n)
      92             : 
      93         583 :                IF (input%l_f) THEN
      94             :                   CALL writeXMLElementFormPoly('forceTotal',(/'atomType',&
      95             :                        'x       ','y       ','z       ',&
      96             :                        'F_x     ','F_y     ','F_z     '/),attributes,&
      97         480 :                         RESHAPE((/8,1,1,1,3,3,3,6,12,12,12,12,12,12/),(/7,2/)))
      98             :                END IF
      99             :             END IF
     100             :          END DO
     101         334 :          IF (input%l_f) CALL closeXMLElement('totalForcesOnRepresentativeAtoms')
     102             : 
     103             :          ! Check convergence of force by comparing force with old_force
     104        2674 :          maxAbsForceDist=MAXVAL(ABS(forcetot - results%force_old))
     105        2674 :          results%force_old(:,:)=forcetot !Store for next iteration
     106        4495 :          results%force=0.0
     107             : 
     108         334 :          l_forceConverged=maxAbsForceDist<input%force_converged
     109         334 :          l_forceConverged=l_forceConverged.and.(results%last_distance.LE.input%mindistance)
     110         334 :          l_forceConverged=l_forceConverged.and.(results%last_distance.GE.0.0)
     111             :          ! In the first iteration results%last_distance is initialized to -1.0.
     112             : 
     113         334 :          IF (.NOT.l_forceConverged) THEN
     114         326 :             WRITE (oUnit,8020) input%force_converged,maxAbsForceDist
     115             : 8020        FORMAT ('No new postions, force convergence required=',f8.5,'; max force distance=',f8.5)
     116             :          END IF
     117             :       END IF
     118             : 
     119             : #ifdef CPP_MPI
     120         668 :       CALL MPI_BCAST(l_forceConverged,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
     121             : #endif
     122             : 
     123         668 :       IF (l_forceConverged.AND.input%l_f.AND.(input%f_level.GE.0).AND.(fmpi%irank==0)) THEN
     124             :          CALL rotate_forces(atoms%ntype,atoms%ntype,atoms%nat,sym%nop,results%tote,&
     125             :                             cell%omtil,atoms%neq,sym%mrot,cell%amat,cell%bmat,&
     126           1 :                             atoms%taual,sym%tau,forcetot,atoms%label)
     127             :       END IF
     128             : 
     129             :     
     130             : 
     131         668 :       IF (l_forceConverged.AND.input%l_f) CALL relaxation(fmpi,input,atoms,cell,sym ,vacuum,forcetot,results%tote)
     132             : 
     133         664 :    END SUBROUTINE force_w
     134             : 
     135             : END MODULE m_forcew