relax_io.F90 7.4 KB
Newer Older
1 2 3 4 5 6 7
!--------------------------------------------------------------------------------
! Copyright (c) 2016 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_relaxio
8 9 10
  !This module handles IO to the relax.xml file
  !The writing is done directly to relax.xml
  !The reading uses the libxml interface to inp.xml. Hence the relax.xml has to be included here.
11 12 13
  USE m_judft
  IMPLICIT NONE
  PRIVATE
14
  PUBLIC :: read_relax,write_relax,apply_displacements,read_displacements
15 16 17 18 19 20 21 22 23 24
CONTAINS
  SUBROUTINE write_relax(positions,forces,energies,displace)
    REAL,INTENT(in):: positions(:,:,:)
    REAL,INTENT(in):: forces(:,:,:)
    REAL,INTENT(in):: energies(:)
    REAL,INTENT(in):: displace(:,:)

    INTEGER :: no_steps,n,ntype,step
    No_steps=SIZE(positions,3)
    ntype=SIZE(positions,2)
25
    IF (ntype.NE.SIZE(forces,2).OR.ntype.NE.SIZE(displace,2).OR.&
26 27 28 29
         no_steps.NE.SIZE(forces,3).OR.no_steps.NE.SIZE(energies))THEN
       CALL judft_error("BUG in relax_io")
    ENDIF
    OPEN(765,file="relax.xml",status="replace")
30 31 32 33
    WRITE(765,*) "<relaxation>"
    !write current set of displacements
    WRITE(765,*) "  <displacements>"
    DO n=1,SIZE(displace,2)
34 35
       WRITE(765,"(a,3(f15.10,1x),a)") &
            '    <displace>',displace(:,n),'</displace>'
36 37
    END DO
    WRITE(765,"(a)") '  </displacements>'
38

39 40 41 42 43
    !Write all known old positions,forces and energies
    WRITE(765,*) "  <relaxation-history>"
    DO step=1,no_steps
       WRITE(765,"(a,f20.10,a)") '    <step energy="',energies(step),'">'
       DO n=1,ntype
44 45
          WRITE(765,"(a,6(f15.10,1x),a)") &
               '      <posforce>',positions(:,n,step),forces(:,n,step),'</posforce>'
46 47 48 49 50 51 52
       END DO
       WRITE(765,"(a)") '    </step>'
    ENDDO
    WRITE(765,*) "  </relaxation-history>"
    WRITE(765,*) "</relaxation>"
    CLOSE(765)
  END SUBROUTINE write_relax
53

54 55 56
  SUBROUTINE read_relax(positions,forces,energies)
    USE m_xmlIntWrapFort 
    USE m_calculator
57 58 59
    REAL,INTENT(INOUT),ALLOCATABLE:: positions(:,:,:)
    REAL,INTENT(INOUT),ALLOCATABLE:: forces(:,:,:)
    REAL,INTENT(INOUT),ALLOCATABLE:: energies(:)
60

61
    REAL,ALLOCATABLE::rtmp(:,:,:)
62 63
    INTEGER:: no_steps
    INTEGER:: ntype,step,n
64
    CHARACTER(len=100):: path,p,str
65
    no_steps=xmlGetNumberOfNodes('/fleurInput/relaxation/relaxation-history/step')
66
    ntype=SIZE(positions,2)
67
    IF (no_steps==0) THEN
68
       IF (.NOT.ALLOCATED(positions)) ALLOCATE(positions(0,0,0),forces(0,0,0),energies(0))
69
       RETURN
70 71 72 73
    END IF
    IF (ALLOCATED(positions)) THEN 
       !Assume that we got already a set of positions, forces, energy and extend that list
       rtmp=positions
74
       DEALLOCATE(positions)
75 76 77
       ALLOCATE(positions(3,ntype,no_steps+SIZE(rtmp,3)))
       positions(:,:,no_steps+1:)=rtmp
       rtmp=forces
78
       DEALLOCATE(forces)
79 80 81
       ALLOCATE(forces(3,ntype,no_steps+SIZE(rtmp,3)))
       forces(:,:,no_steps+1:)=rtmp
       rtmp(1,1,:)=energies
82
       DEALLOCATE(energies)
83 84
       ALLOCATE(energies(no_steps+SIZE(rtmp,3)))
       energies(no_steps+1:)=rtmp(1,1,:)
85
    ELSE
86 87 88
       ALLOCATE(positions(3,ntype,no_steps))
       ALLOCATE(forces(3,ntype,no_steps))
       ALLOCATE(energies(no_steps))
89 90 91
    END IF
    DO step=1,no_steps
       WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/relaxation-history/step[',step,']'
92
       energies(step)=evaluateFirstOnly(xmlGetAttributeValue(TRIM(path)//"/@energy"))
93
       DO n=1,ntype
94
          WRITE(p,"(a,a,i0,a)") TRIM(path),"/posforce[",n,"]"
95 96
          str=xmlGetAttributeValue(p)
          positions(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
97
          Forces(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
98 99 100 101
       ENDDO
    END DO
  END SUBROUTINE read_relax

102

103
  SUBROUTINE read_displacements(atoms,disp)
104 105 106
    USE m_xmlIntWrapFort 
    USE m_calculator
    USE m_types
107 108 109
    TYPE(t_atoms),INTENT(in)::atoms
    REAL,INTENT(out)::disp(:,:)
    CHARACTER(len=50):: path,str
110 111 112 113 114 115 116 117 118 119 120 121 122
    INTEGER :: n
    disp=0.0
    IF (xmlGetNumberOfNodes('/fleurInput/relaxation/displacements')==0) RETURN
    !read displacements and apply to positions
    IF (atoms%ntype.NE.xmlGetNumberOfNodes('/fleurInput/relaxation/displacements/displace')) CALL judft_error("Wrong number of displacements in relaxation")
    DO n=1,atoms%ntype
       WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/displacements/displace[',n,']'
       str=xmlGetAttributeValue(path)
       disp(:,n)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
    END DO
  END SUBROUTINE read_displacements

  SUBROUTINE apply_displacements(cell,input,vacuum,oneD,sym,noco,atoms)
123 124
    USE m_types
    USE m_chkmt
125
    USE m_mapatom
126 127 128 129 130 131 132
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_vacuum),INTENT(IN)  :: vacuum
    TYPE(t_cell),INTENT(IN)    :: cell
    TYPE(t_oneD),INTENT(IN)    :: oneD
    TYPE(t_sym),INTENT(INOUT)  :: sym
    TYPE(t_noco),INTENT(IN)    :: noco

133
    TYPE(t_atoms),INTENT(INOUT):: atoms
134 135


136 137
    INTEGER:: n,indx(2)
    REAL   :: disp(3,atoms%ntype),disp_all(3,atoms%nat),taual0(3,atoms%nat),overlap(0:atoms%ntype,atoms%ntype)
138

139 140 141 142 143
    CALL read_displacements(atoms,disp)
    IF (ALL(ABS(disp)<1E-8)) RETURN
    !Fist make sure original MT spheres do not overlap
    CALL chkmt(atoms,input,vacuum,cell,oneD,.TRUE.,overlap=overlap)
    IF(ANY(overlap>0.0)) CALL judft_error("Overlapping MT-spheres in relaxation before even applying any displacement",hint="You messed up your setup")
144

145
    taual0=atoms%taual !Store original positions
146

147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
    !Now check for overlapping mt-spheres
    overlap=1.0
    DO WHILE(ANY(overlap>1E-10))
       atoms%taual=taual0  
       CALL rotate_to_all_sites(disp,atoms,cell,sym,disp_all)
       atoms%taual=taual0+disp_all
       atoms%pos=MATMUL(cell%amat,atoms%taual)
       CALL chkmt(atoms,input,vacuum,cell,oneD,.TRUE.,overlap=overlap)
       IF (ANY(overlap>0.0)) THEN
          IF (ANY(overlap(0,:)>1E-10)) CALL judft_error("Atom spills out into vacuum during relaxation")
          indx=MAXLOC(overlap(1:,:)) !the two indices with the most overlap
          !Try only 90% of displacement
          disp(:,indx(1))=disp(:,indx(1))*0.9 
          disp(:,indx(2))=disp(:,indx(2))*0.9
          WRITE(*,*) "Attention: Overlapping MT-spheres. Reduced displacement by 10%"
          WRITE(*,*) indx,overlap(indx(1),indx(2))
       END IF
    END DO
165

166
    CALL mapatom(sym,atoms,cell,input,noco)
167

168 169 170 171 172 173 174 175 176 177 178 179 180 181
    WRITE(6,*) "Atomic positions including displacements:"
    DO n=1,atoms%nat
       WRITE(6,*) n,atoms%taual(:,n),atoms%pos(:,n)
    ENDDO

  END SUBROUTINE apply_displacements

  SUBROUTINE rotate_to_all_sites(disp,atoms,cell,sym,disp_all)
    USE m_types
    REAL,INTENT(in)          :: disp(:,:)
    TYPE(t_atoms),INTENT(in) :: atoms
    TYPE(t_cell),INTENT(IN)  :: cell
    TYPE(t_sym),INTENT(IN)   :: sym
    REAL,INTENT(out)         :: disp_all(:,:)
182

183 184
    INTEGER:: n,na,jop
    REAL   :: tau0(3),tau0_rot(3),tau_rot(3)
185 186


187 188 189 190 191 192 193
    DO n=1,atoms%ntype
       tau0=atoms%taual(:,n)
       DO na=SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
          jop = sym%invtab(atoms%ngopr(na))
          tau0_rot=MATMUL(1.*sym%mrot(:,:,jop),tau0)+sym%tau(:,jop) !translation will cancel, included for clarity
          tau_rot=MATMUL(1.*sym%mrot(:,:,jop),tau0+disp(:,n))+sym%tau(:,jop)
          disp_all(:,na)=tau_rot-tau0_rot
194
       END DO
195 196
    END DO
  END SUBROUTINE rotate_to_all_sites
197
END MODULE m_relaxio