force_w.f90 4.36 KB
Newer Older
1 2 3 4 5 6
      MODULE m_forcew
! ************************************************************
! Printing force components
! ************************************************************
      CONTAINS
      SUBROUTINE force_w(&
7
     &                   input,atoms,sym,results,cell,oneD,vacuum)
8 9 10
      USE m_geo
      USE m_relax
      USE m_types
11
      USE m_xmlOutput
12 13 14 15 16 17 18 19
      IMPLICIT NONE

      TYPE(t_results),INTENT(IN)   :: 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
20
      TYPE(t_vacuum),INTENT(IN)    :: vacuum
21 22 23 24 25 26 27 28 29
!     ..
!     .. Local Scalars ..
      REAL,PARAMETER:: zero=0.0
      REAL sum
      INTEGER i,jsp,n,nat1,ierr
      REAL eps_force
      LOGICAL :: l_new
!     ..
!     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
30
      REAL forcetot(3,atoms%ntype)
31
      CHARACTER(LEN=20) :: attributes(7)
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
!
!     write spin-dependent forces
!
      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)
               WRITE (16,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)
      WRITE (16,8005)
 8005 FORMAT (/,' ***** TOTAL FORCES ON ATOMS ***** ',/)
58
      IF (input%l_f) CALL openXMLElement('totalForcesOnRepresentativeAtoms',(/'units'/),(/'Htr/bohr'/))
59 60 61 62
      nat1 = 1
      DO n = 1,atoms%ntype
         IF (atoms%l_geo(n)) THEN
            DO i = 1,3
63
               forcetot(i,n) = zero
64
            END DO
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
            DO jsp = 1,input%jspins
               DO i = 1,3
                  forcetot(i,n) = forcetot(i,n) + results%force(i,n,jsp)
               END DO
            END DO

            WRITE (6,FMT=8010) n, (atoms%pos(i,nat1),i=1,3),&
     &        (forcetot(i,n),i=1,3)
            WRITE (16,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)

            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
92 93
         nat1 = nat1 + atoms%neq(n)
      END DO
94
      IF (input%l_f) CALL closeXMLElement('totalForcesOnRepresentativeAtoms')
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

      sum=0.0
      DO n = 1,atoms%ntype
        IF (atoms%l_geo(n)) THEN
          DO i = 1,3
            sum = max(sum,(forcetot(i,n) - results%force_old(i,n))**2)
          ENDDO
        ENDIF
      ENDDO 
      sum=sqrt(sum)
!-roa
      eps_force=0.00001
      open(88,file='eps_force',form='formatted',status='old',err=188)
      read(88,'(f20.8)') eps_force
      close (88)
  188 continue
!-roa
 
      WRITE (6,8020) eps_force,sum
 8020 FORMAT ('eps_force=',f8.5,'max=',f8.5)

      INQUIRE(file ="relax_inp",exist= l_new)
      IF (l_new) THEN
        CALL relax(input%film,atoms%pos,atoms%neq,sym%mrot,sym%tau,cell%amat,cell%bmat,atoms%ngopr,sym%invtab&
     &        ,forcetot)
      ELSE

         IF ((sum<eps_force).AND.input%l_f) THEN
            CALL geo(&
     &           atoms,sym,cell,&
125
     &           oneD,vacuum,input,results%tote,forcetot)
126 127 128 129 130
         END IF
      ENDIF

      END SUBROUTINE force_w
      END MODULE m_forcew