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
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 55 56
  SUBROUTINE read_relax(positions,forces,energies)
    USE m_xmlIntWrapFort 
    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
    no_steps=xmlGetNumberOfNodes('/fleurInput/relaxation/relaxation-history/step')
Daniel Wortmann's avatar
Daniel Wortmann committed
66
    ntype=SIZE(positions,2)
67
    IF (no_steps==0) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
68
       IF (.NOT.ALLOCATED(positions)) ALLOCATE(positions(0,0,0),forces(0,0,0),energies(0))
69
       RETURN
Daniel Wortmann's avatar
Daniel Wortmann committed
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)
Daniel Wortmann's avatar
Daniel Wortmann committed
75 76 77
       ALLOCATE(positions(3,ntype,no_steps+SIZE(rtmp,3)))
       positions(:,:,no_steps+1:)=rtmp
       rtmp=forces
78
       DEALLOCATE(forces)
Daniel Wortmann's avatar
Daniel Wortmann committed
79 80 81
       ALLOCATE(forces(3,ntype,no_steps+SIZE(rtmp,3)))
       forces(:,:,no_steps+1:)=rtmp
       rtmp(1,1,:)=energies
82
       DEALLOCATE(energies)
Daniel Wortmann's avatar
Daniel Wortmann committed
83 84
       ALLOCATE(energies(no_steps+SIZE(rtmp,3)))
       energies(no_steps+1:)=rtmp(1,1,:)
85
    ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
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,']'
Daniel Wortmann's avatar
Daniel Wortmann committed
92
       energies(step)=evaluateFirstOnly(xmlGetAttributeValue(TRIM(path)//"/@energy"))
93
       DO n=1,ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
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)/)
Daniel Wortmann's avatar
Daniel Wortmann committed
97
          Forces(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
98 99 100 101
       ENDDO
    END DO
  END SUBROUTINE read_relax

102

Daniel Wortmann's avatar
Daniel Wortmann committed
103
  SUBROUTINE read_displacements(atoms,disp)
104 105 106
    USE m_xmlIntWrapFort 
    USE m_calculator
    USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
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)
Daniel Wortmann's avatar
Daniel Wortmann committed
123 124
    USE m_types
    USE m_chkmt
125
    USE m_mapatom
Daniel Wortmann's avatar
Daniel Wortmann committed
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


Daniel Wortmann's avatar
Daniel Wortmann committed
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

Daniel Wortmann's avatar
Daniel Wortmann committed
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

Daniel Wortmann's avatar
Daniel Wortmann committed
145
    taual0=atoms%taual !Store original positions
146

Daniel Wortmann's avatar
Daniel Wortmann committed
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

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

Daniel Wortmann's avatar
Daniel Wortmann committed
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

Daniel Wortmann's avatar
Daniel Wortmann committed
183 184
    INTEGER:: n,na,jop
    REAL   :: tau0(3),tau0_rot(3),tau_rot(3)
185 186


Daniel Wortmann's avatar
Daniel Wortmann committed
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
Daniel Wortmann's avatar
Daniel Wortmann committed
195 196
    END DO
  END SUBROUTINE rotate_to_all_sites
197
END MODULE m_relaxio