Line data    Source code

       1             : MODULE m_rotate_forces
       2             :    ! This routine writes a file similar to forces.dat, but containing
       3             :    ! all atoms in a sequence corresponding to their appearance in the
       4             :    ! input file for the input-file generator. forces.dat only
       5             :    ! contains forces for a representant of the symmetry equivalent
       6             :    ! atoms. This routine also gives files FORCES and POSCAR to be
       7             :    ! used with the phon package by Dario Alfe (for a single
       8             :    ! displacement),
       9             :    ! Reference
      10             :    ! Klueppelberg, Jun 2012
      11             : 
      12             :    ! Readded to only construct POSCAR and FORCES files in Dec 2020, Neukirchen
      13             : 
      14             :    ! Modified to construct a file for use with phonopy instead of PHON.
      15             :    ! Neukirchen, Dec 2020 
      16             : 
      17             : CONTAINS
      18           1 :    SUBROUTINE rotate_forces(ntypd,ntype,natd,nop,tote,omtil,neq,mrot,amat,bmat,taual,tau,force,label)
      19             :       USE m_constants
      20             :       USE m_juDFT_string
      21             : 
      22             :       IMPLICIT NONE
      23             : 
      24             :       INTEGER, INTENT (IN) :: ntypd,ntype,natd,nop
      25             :       REAL, INTENT (IN) :: tote,omtil
      26             : 
      27             :       INTEGER, INTENT (IN) :: neq(ntypd),mrot(3,3,nop)
      28             :       REAL, INTENT (IN) :: amat(3,3),bmat(3,3),taual(3,natd),tau(3,nop)
      29             :       REAL, INTENT (IN) :: force(3,ntypd)
      30             :       CHARACTER(LEN=20) :: label(natd)
      31             : 
      32           1 :       INTEGER :: iatom,iatom0,itype,ieq,ratom,isym,sym,i,true_atom,atom_map(natd)
      33             : 
      34           1 :       REAL :: pos(3,natd),rtaual(3),forcerot(3,3),forceval(3,natd)
      35             :       CHARACTER(len=20) :: string
      36             :       LOGICAL :: l_PHON
      37             : 
      38           1 :       l_PHON = .FALSE.
      39           1 :       OPEN (79,file='FORCES')
      40             :       IF (l_PHON) THEN
      41             :          OPEN (80,file='POSCAR')
      42             :       END IF
      43           1 :       OPEN(81,file='FORCES_SORT')
      44             : 
      45           1 :       WRITE (79,'(i1)') 1
      46           1 :       WRITE (79,'(i1,1x,a)') 1,'#'
      47           1 :       WRITE (81,'(i1)') 1
      48           1 :       WRITE (81,'(i1,1x,a)') 1,'#'
      49             :       IF (l_PHON) THEN
      50             :          WRITE (80,'(a)') "#insert lattice name"
      51             :          WRITE (80,'(f20.10)') -omtil*bohr_to_angstrom_const**3
      52             :          WRITE (80,FMT=800) amat(1,1:3)*bohr_to_angstrom_const
      53             :          WRITE (80,FMT=800) amat(2,1:3)*bohr_to_angstrom_const
      54             :          WRITE (80,FMT=800) amat(3,1:3)*bohr_to_angstrom_const
      55             :          WRITE (string,'(i3)') ntype
      56             :          WRITE (80,'('//string//'(i2,1x))') (neq(itype),itype=1,ntype)
      57             :          WRITE (80,'(a)') 'Direct'
      58             :       END IF
      59             : 
      60           1 :       iatom  = 0
      61           1 :       iatom0 = 1
      62           3 :       DO itype = 1,ntype
      63           5 :          DO ieq = 1,neq(itype)
      64           3 :             iatom = iatom+1
      65           3 :             true_atom = str2int(TRIM(label(iatom)))
      66           3 :             atom_map(true_atom) = iatom
      67          48 :             pos(:,iatom) = matmul(amat,taual(:,iatom))
      68             : 
      69             :             ratom = 0
      70           4 :             sym = 0
      71             : 
      72           4 :             DO isym = 1,nop
      73         112 :                rtaual(:)=matmul(mrot(:,:,isym),taual(:,iatom0))+tau(:,isym)
      74          13 :                IF(all(abs(modulo(rtaual-taual(:,iatom)+0.5,1d0)-0.5).lt.1d-10)) THEN
      75           3 :                   sym = isym
      76           3 :                   EXIT
      77             :                END IF
      78             :             END DO
      79             : 
      80           3 :             IF (sym.eq.0) THEN
      81           0 :                sym = 1
      82             :             END IF
      83             : 
      84         309 :             forcerot = matmul(amat,matmul(mrot(:,:,sym),bmat/tpi_const))
      85          48 :             forceval(:,iatom) = matmul(forcerot,force(:,itype))
      86             : 
      87           2 :             IF (l_PHON) THEN
      88             :                WRITE (79,FMT=790) forceval(1:3,iatom)*hartree_to_ev_const/bohr_to_angstrom_const
      89             :                WRITE (80,FMT=800) taual(1:3,iatom)
      90             :             ELSE
      91           3 :                WRITE (79,*) forceval(1:3,iatom), 'force'
      92             :             END IF
      93             :          END DO
      94           3 :          iatom0 = iatom0 + neq(itype)
      95             :       END DO
      96             : 
      97           4 :       DO iatom = 1, natd
      98           4 :          WRITE (81,*) forceval(1:3,atom_map(iatom)), 'force'
      99             :       END DO
     100             : 
     101             : 790   FORMAT (1x,3(1x,f20.10))
     102             : 800   FORMAT (1x,3(1x,f20.10))
     103             : 810   FORMAT (1x,a4,43x,3f14.9)
     104             : 
     105           1 :       CLOSE(79)
     106             :       IF (l_PHON) THEN
     107             :          CLOSE(80)
     108             :       END IF
     109           1 :       CLOSE(81)
     110             : 
     111           1 :    END SUBROUTINE rotate_forces
     112             : 
     113             : END MODULE m_rotate_forces