relax_io.F90 7.45 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
Daniel Wortmann's avatar
Daniel Wortmann committed
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.&
Daniel Wortmann's avatar
Daniel Wortmann committed
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)
Daniel Wortmann's avatar
Daniel Wortmann committed
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
Daniel Wortmann's avatar
Daniel Wortmann committed
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
  SUBROUTINE read_relax(positions,forces,energies)
55
    USE m_types_xml 
56
    USE m_calculator
Daniel Wortmann's avatar
Daniel Wortmann committed
57 58 59
    REAL,INTENT(INOUT),ALLOCATABLE:: positions(:,:,:)
    REAL,INTENT(INOUT),ALLOCATABLE:: forces(:,:,:)
    REAL,INTENT(INOUT),ALLOCATABLE:: energies(:)
60

Daniel Wortmann's avatar
Daniel Wortmann committed
61
    REAL,ALLOCATABLE::rtmp(:,:,:)
62 63
    INTEGER:: no_steps
    INTEGER:: ntype,step,n
Daniel Wortmann's avatar
Daniel Wortmann committed
64
    CHARACTER(len=100):: path,p,str
65 66 67 68

    TYPE(t_xml)::xml

    no_steps=xml%GetNumberOfNodes('/fleurInput/relaxation/relaxation-history/step')
Daniel Wortmann's avatar
Daniel Wortmann committed
69
    ntype=SIZE(positions,2)
70
    IF (no_steps==0) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
71
       IF (.NOT.ALLOCATED(positions)) ALLOCATE(positions(0,0,0),forces(0,0,0),energies(0))
72
       RETURN
Daniel Wortmann's avatar
Daniel Wortmann committed
73 74 75 76
    END IF
    IF (ALLOCATED(positions)) THEN 
       !Assume that we got already a set of positions, forces, energy and extend that list
       rtmp=positions
77
       DEALLOCATE(positions)
Daniel Wortmann's avatar
Daniel Wortmann committed
78 79 80
       ALLOCATE(positions(3,ntype,no_steps+SIZE(rtmp,3)))
       positions(:,:,no_steps+1:)=rtmp
       rtmp=forces
81
       DEALLOCATE(forces)
Daniel Wortmann's avatar
Daniel Wortmann committed
82 83 84
       ALLOCATE(forces(3,ntype,no_steps+SIZE(rtmp,3)))
       forces(:,:,no_steps+1:)=rtmp
       rtmp(1,1,:)=energies
85
       DEALLOCATE(energies)
Daniel Wortmann's avatar
Daniel Wortmann committed
86 87
       ALLOCATE(energies(no_steps+SIZE(rtmp,3)))
       energies(no_steps+1:)=rtmp(1,1,:)
88
    ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
89 90 91
       ALLOCATE(positions(3,ntype,no_steps))
       ALLOCATE(forces(3,ntype,no_steps))
       ALLOCATE(energies(no_steps))
92 93 94
    END IF
    DO step=1,no_steps
       WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/relaxation-history/step[',step,']'
95
       energies(step)=evaluateFirstOnly(xml%GetAttributeValue(TRIM(path)//"/@energy"))
96
       DO n=1,ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
97
          WRITE(p,"(a,a,i0,a)") TRIM(path),"/posforce[",n,"]"
98
          str=xml%GetAttributeValue(p)
99
          positions(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
Daniel Wortmann's avatar
Daniel Wortmann committed
100
          Forces(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
101 102 103 104
       ENDDO
    END DO
  END SUBROUTINE read_relax

105

Daniel Wortmann's avatar
Daniel Wortmann committed
106
  SUBROUTINE read_displacements(atoms,disp)
107
    USE m_types_xml 
108 109
    USE m_calculator
    USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
110 111 112
    TYPE(t_atoms),INTENT(in)::atoms
    REAL,INTENT(out)::disp(:,:)
    CHARACTER(len=50):: path,str
113
    INTEGER :: n
114 115 116

    TYPE(t_xml)::xml

117
    disp=0.0
118
    IF (xml%GetNumberOfNodes('/fleurInput/relaxation/displacements')==0) RETURN
119
    !read displacements and apply to positions
120
    IF (atoms%ntype.NE.xml%GetNumberOfNodes('/fleurInput/relaxation/displacements/displace')) CALL judft_error("Wrong number of displacements in relaxation")
121 122
    DO n=1,atoms%ntype
       WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/displacements/displace[',n,']'
123
       str=xml%GetAttributeValue(path)
124 125 126 127 128
       disp(:,n)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
    END DO
  END SUBROUTINE read_displacements

  SUBROUTINE apply_displacements(cell,input,vacuum,oneD,sym,noco,atoms)
Daniel Wortmann's avatar
Daniel Wortmann committed
129 130
    USE m_types
    USE m_chkmt
131
    USE m_mapatom
Daniel Wortmann's avatar
Daniel Wortmann committed
132 133 134 135 136 137 138
    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

139
    TYPE(t_atoms),INTENT(INOUT):: atoms
140 141


Daniel Wortmann's avatar
Daniel Wortmann committed
142 143
    INTEGER:: n,indx(2)
    REAL   :: disp(3,atoms%ntype),disp_all(3,atoms%nat),taual0(3,atoms%nat),overlap(0:atoms%ntype,atoms%ntype)
144

Daniel Wortmann's avatar
Daniel Wortmann committed
145 146 147 148 149
    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")
150

Daniel Wortmann's avatar
Daniel Wortmann committed
151
    taual0=atoms%taual !Store original positions
152

Daniel Wortmann's avatar
Daniel Wortmann committed
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
    !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
171

Daniel Wortmann's avatar
Daniel Wortmann committed
172
    CALL mapatom(sym,atoms,cell,input,noco)
173

Daniel Wortmann's avatar
Daniel Wortmann committed
174 175 176 177 178 179 180 181 182 183 184 185 186 187
    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(:,:)
188

Daniel Wortmann's avatar
Daniel Wortmann committed
189 190
    INTEGER:: n,na,jop
    REAL   :: tau0(3),tau0_rot(3),tau_rot(3)
191 192


Daniel Wortmann's avatar
Daniel Wortmann committed
193 194 195 196 197 198 199
    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
200
       END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
201 202
    END DO
  END SUBROUTINE rotate_to_all_sites
203
END MODULE m_relaxio