Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_relaxio
8 : !This module handles IO to the relax.xml file
9 : !The writing is done directly to relax.xml
10 : !The reading uses the libxml interface to inp.xml. Hence the relax.xml has to be included here.
11 : USE m_judft
12 : USE m_constants
13 : IMPLICIT NONE
14 : PRIVATE
15 : PUBLIC :: read_relax,write_relax,apply_displacements,read_displacements,rotate_to_all_sites
16 : CONTAINS
17 2 : SUBROUTINE write_relax(positions,forces,energies,displace)
18 : REAL,INTENT(in):: positions(:,:,:)
19 : REAL,INTENT(in):: forces(:,:,:)
20 : REAL,INTENT(in):: energies(:)
21 : REAL,INTENT(in):: displace(:,:)
22 :
23 : INTEGER :: no_steps,n,ntype,step
24 :
25 : CHARACTER(len=100):: path,p,str,filename_add
26 :
27 2 : filename_add = ""
28 2 : IF (judft_was_argument("-add_name")) filename_add = TRIM(judft_string_for_argument("-add_name"))//"_"
29 :
30 2 : No_steps=SIZE(positions,3)
31 2 : ntype=SIZE(positions,2)
32 : IF (ntype.NE.SIZE(forces,2).OR.ntype.NE.SIZE(displace,2).OR.&
33 2 : no_steps.NE.SIZE(forces,3).OR.no_steps.NE.SIZE(energies))THEN
34 0 : CALL judft_error("BUG in relax_io")
35 : ENDIF
36 2 : OPEN(765,file=TRIM(filename_add)//"relax.xml",status="replace")
37 2 : WRITE(765,*) "<!-- Attention, absolute coordinates used here -->"
38 2 : WRITE(765,*) "<relaxation>"
39 : !write current set of displacements
40 2 : WRITE(765,*) " <displacements>"
41 6 : DO n=1,SIZE(displace,2)
42 : WRITE(765,"(a,3(f15.10,1x),a)") &
43 6 : ' <displace>',displace(:,n),'</displace>'
44 : END DO
45 2 : WRITE(765,"(a)") ' </displacements>'
46 :
47 : !Write all known old positions,forces and energies
48 2 : WRITE(765,*) " <relaxation-history>"
49 4 : DO step=1,no_steps
50 2 : WRITE(765,"(a,f20.10,a)") ' <step energy="',energies(step),'">'
51 6 : DO n=1,ntype
52 : WRITE(765,"(a,6(f15.10,1x),a)") &
53 6 : ' <posforce>',positions(:,n,step),forces(:,n,step),'</posforce>'
54 : END DO
55 4 : WRITE(765,"(a)") ' </step>'
56 : ENDDO
57 2 : WRITE(765,*) " </relaxation-history>"
58 2 : WRITE(765,*) "</relaxation>"
59 2 : CLOSE(765)
60 2 : END SUBROUTINE write_relax
61 :
62 2 : SUBROUTINE read_relax(positions,forces,energies)
63 : USE m_types_xml
64 : USE m_calculator
65 : REAL,INTENT(INOUT),ALLOCATABLE:: positions(:,:,:)
66 : REAL,INTENT(INOUT),ALLOCATABLE:: forces(:,:,:)
67 : REAL,INTENT(INOUT),ALLOCATABLE:: energies(:)
68 :
69 2 : REAL,ALLOCATABLE::rtmp(:,:,:)
70 : INTEGER:: no_steps
71 : INTEGER:: ntype,step,n
72 : CHARACTER(len=100):: path,p,str,filename_add
73 :
74 : TYPE(t_xml)::xml
75 2 : filename_add = ""
76 2 : IF (judft_was_argument("-add_name")) filename_add = TRIM(judft_string_for_argument("-add_name"))//"_"
77 2 : CALL xml%init(filename_add)
78 2 : no_steps=xml%GetNumberOfNodes('/fleurInput/relaxation/relaxation-history/step')
79 2 : ntype=SIZE(positions,2)
80 2 : IF (no_steps==0) THEN
81 2 : IF (.NOT.ALLOCATED(positions)) ALLOCATE(positions(0,0,0),forces(0,0,0),energies(0))
82 : RETURN
83 : END IF
84 0 : IF (ALLOCATED(positions)) THEN
85 : !Assume that we got already a set of positions, forces, energy and extend that list
86 0 : rtmp=positions
87 0 : DEALLOCATE(positions)
88 0 : ALLOCATE(positions(3,ntype,no_steps+SIZE(rtmp,3)))
89 0 : positions(:,:,no_steps+1:)=rtmp
90 0 : rtmp=forces
91 0 : DEALLOCATE(forces)
92 0 : ALLOCATE(forces(3,ntype,no_steps+SIZE(rtmp,3)))
93 0 : forces(:,:,no_steps+1:)=rtmp
94 0 : rtmp(1,1,:)=energies
95 0 : DEALLOCATE(energies)
96 0 : ALLOCATE(energies(no_steps+SIZE(rtmp,3)))
97 0 : energies(no_steps+1:)=rtmp(1,1,:)
98 : ELSE
99 0 : ALLOCATE(positions(3,ntype,no_steps))
100 0 : ALLOCATE(forces(3,ntype,no_steps))
101 0 : ALLOCATE(energies(no_steps))
102 : END IF
103 0 : DO step=1,no_steps
104 0 : WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/relaxation-history/step[',step,']'
105 0 : energies(step)=evaluateFirstOnly(xml%GetAttributeValue(TRIM(path)//"/@energy"))
106 0 : DO n=1,ntype
107 0 : WRITE(p,"(a,a,i0,a)") TRIM(path),"/posforce[",n,"]"
108 0 : str=xml%GetAttributeValue(p)
109 0 : positions(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
110 0 : Forces(:,n,step)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
111 : ENDDO
112 : END DO
113 0 : call xml%FreeResources()
114 2 : END SUBROUTINE read_relax
115 :
116 :
117 80 : SUBROUTINE read_displacements(atoms,disp)
118 : USE m_types_xml
119 : USE m_calculator
120 : USE m_types
121 : TYPE(t_atoms),INTENT(in)::atoms
122 : REAL,INTENT(out)::disp(:,:)
123 : CHARACTER(len=50):: path,str
124 : INTEGER :: n
125 :
126 : TYPE(t_xml)::xml
127 :
128 628 : disp=0.0
129 80 : IF (xml%GetNumberOfNodes('/fleurInput/relaxation/displacements')==0) RETURN
130 : !read displacements and apply to positions
131 0 : IF (atoms%ntype.NE.xml%GetNumberOfNodes('/fleurInput/relaxation/displacements/displace')) CALL judft_error("Wrong number of displacements in relaxation")
132 0 : DO n=1,atoms%ntype
133 0 : WRITE(path,"(a,i0,a)") '/fleurInput/relaxation/displacements/displace[',n,']'
134 0 : str=xml%GetAttributeValue(path)
135 0 : disp(:,n)=(/evaluateFirst(str),evaluateFirst(str),evaluateFirst(str)/)
136 : END DO
137 80 : END SUBROUTINE read_displacements
138 :
139 80 : SUBROUTINE apply_displacements(cell,input,vacuum ,sym,noco,atoms,gfinp)
140 : USE m_types
141 : USE m_chkmt
142 : USE m_constants
143 : USE m_mapatom
144 : TYPE(t_input),INTENT(IN) :: input
145 : TYPE(t_vacuum),INTENT(IN) :: vacuum
146 : TYPE(t_cell),INTENT(IN) :: cell
147 :
148 : TYPE(t_sym),INTENT(INOUT) :: sym
149 : TYPE(t_noco),INTENT(IN) :: noco
150 : type(t_gfinp), intent(in) :: gfinp
151 80 : character(len=:), allocatable :: error_output
152 :
153 : TYPE(t_atoms),INTENT(INOUT):: atoms
154 :
155 :
156 : INTEGER:: n,indx(2)
157 80 : REAL :: disp(3,atoms%ntype),disp_all(3,atoms%nat),taual0(3,atoms%nat),overlap(0:atoms%ntype,atoms%ntype)
158 :
159 80 : CALL read_displacements(atoms,disp)
160 : !change displacements to internal coordinates
161 3197 : disp=MATMUL(cell%bmat,disp)/tpi_const
162 628 : IF (ALL(ABS(disp)<1E-8)) RETURN
163 : !Fist make sure original MT spheres do not overlap
164 0 : CALL chkmt(atoms,input,vacuum,cell ,.TRUE.,overlap=overlap)
165 0 : IF(ANY(overlap>0.0)) CALL judft_error("Overlapping MT-spheres in relaxation before even applying any displacement",hint="You messed up your setup")
166 :
167 0 : taual0=atoms%taual !Store original positions
168 :
169 : !Now check for overlapping mt-spheres
170 0 : overlap=1.0
171 0 : DO WHILE(ANY(overlap>1E-10))
172 0 : atoms%taual=taual0
173 0 : CALL rotate_to_all_sites(disp,atoms,cell,sym,disp_all)
174 0 : atoms%taual=taual0+disp_all
175 0 : atoms%pos=MATMUL(cell%amat,atoms%taual)
176 0 : CALL chkmt(atoms,input,vacuum,cell ,.TRUE.,overlap=overlap)
177 0 : IF (ANY(overlap>0.0)) THEN
178 0 : IF (ANY(overlap(0,:)>1E-10)) CALL judft_error("Atom spills out into vacuum during relaxation")
179 0 : indx=MAXLOC(overlap(1:,:)) !the two indices with the most overlap
180 : !Try only 90% of displacement
181 0 : disp(:,indx(1))=disp(:,indx(1))*0.9
182 0 : disp(:,indx(2))=disp(:,indx(2))*0.9
183 0 : WRITE(*,*) "Attention: Overlapping MT-spheres. Reduced displacement by 10%"
184 0 : WRITE(*,*) indx,overlap(indx(1),indx(2))
185 0 : WRITE(oUnit,'(a,2(i0,1x),f12.8)') "Attention, overlapping MT-spheres: ",indx,overlap(indx(1),indx(2))
186 : ! WRITE(error_output, '(3a,f12.8,a)') "Overlapping MT-spheres during relaxation: ", atoms%label(atoms%firstAtom(indx(1))),&
187 : ! &atoms%label(atoms%firstAtom(indx(2))), overlap(indx(1),indx(2)),&
188 : ! &NEW_LINE('A')//"Treat as an error: writing rescaled displacements to relax.xml is not implemented"
189 : error_output = "Overlapping MT-spheres during relaxation: " // NEW_LINE('A') // &
190 : strip(atoms%label(atoms%firstAtom(indx(1)))) // " " // &
191 : strip(atoms%label(atoms%firstAtom(indx(2)))) // &
192 : " olap: "//float2str(overlap(indx(1),indx(2)))//NEW_LINE('A')//&
193 0 : "Treat as an error: writing rescaled displacements to relax.xml is not implemented"
194 0 : CALL judft_error(error_output)
195 : END IF
196 : END DO
197 :
198 0 : CALL mapatom(sym,atoms,cell,input,noco,gfinp)
199 :
200 0 : WRITE(oUnit,*) "Atomic positions including displacements:"
201 0 : DO n=1,atoms%nat
202 0 : WRITE(oUnit,"(i4,6(1x,f12.5))") n,atoms%taual(:,n),atoms%pos(:,n)
203 : ENDDO
204 :
205 0 : END SUBROUTINE apply_displacements
206 :
207 4 : SUBROUTINE rotate_to_all_sites(disp,atoms,cell,sym,disp_all)
208 : USE m_types
209 : REAL,INTENT(in) :: disp(:,:)
210 : TYPE(t_atoms),INTENT(in) :: atoms
211 : TYPE(t_cell),INTENT(IN) :: cell
212 : TYPE(t_sym),INTENT(IN) :: sym
213 : REAL,INTENT(out) :: disp_all(:,:)
214 :
215 : INTEGER:: iType,iAtom,jop, startAtom
216 : REAL :: tau0(3),tau0_rot(3),tau_rot(3)
217 :
218 :
219 12 : DO iType = 1, atoms%ntype
220 8 : startAtom = atoms%firstAtom(iType)
221 32 : tau0=atoms%taual(:,startAtom)
222 22 : DO iAtom = startAtom, startAtom + atoms%neq(iType) - 1
223 10 : jop = sym%invtab(sym%ngopr(iAtom))
224 280 : tau0_rot=MATMUL(1.*sym%mrot(:,:,jop),tau0)+sym%tau(:,jop) !translation will cancel, included for clarity
225 320 : tau_rot=MATMUL(1.*sym%mrot(:,:,jop),tau0+disp(:,iType))+sym%tau(:,jop)
226 48 : disp_all(:,iAtom)=tau_rot-tau0_rot
227 : END DO
228 : END DO
229 4 : END SUBROUTINE rotate_to_all_sites
230 100 : END MODULE m_relaxio
|
