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_relaxation
8 : USE m_judft
9 :
10 : #ifdef CPP_MPI
11 : USE mpi
12 : #endif
13 :
14 : IMPLICIT NONE
15 :
16 : PRIVATE
17 : PUBLIC relaxation !This is the interface. Below there are internal subroutines for bfgs, simple mixing, CG ...
18 :
19 : CONTAINS
20 4 : SUBROUTINE relaxation(fmpi,input,atoms,cell,sym ,vacuum,force_new,energies_new)
21 : ! This routine uses the current force,energies and atomic positions to
22 : ! generate a displacement in a relaxation step.
23 : ! The history is taken into account by read_relax from m_relaxio
24 : ! After generating new positions the code stops
25 :
26 : USE m_types
27 : USE m_constants
28 : USE m_relaxio
29 : USE m_mixing_history
30 : USE m_chkmt
31 : USE m_types_xml
32 : USE m_xsf_io
33 :
34 : TYPE(t_mpi), INTENT(IN) :: fmpi
35 : TYPE(t_input), INTENT(IN) :: input
36 : TYPE(t_atoms), INTENT(IN) :: atoms
37 : TYPE(t_sym), INTENT(IN) :: sym
38 :
39 : TYPE(t_vacuum), INTENT(IN) :: vacuum
40 : TYPE(t_cell), INTENT(IN) :: cell
41 : REAL, INTENT(IN) :: force_new(:,:), energies_new !data for this iteration
42 :
43 4 : REAL, ALLOCATABLE :: pos(:,:,:), force(:,:,:), energies(:)
44 4 : REAL, ALLOCATABLE :: displace(:,:), old_displace(:,:), tempDisplace(:,:)
45 4 : REAL, ALLOCATABLE :: totalDisplace(:,:)
46 4 : REAL :: dispAll(3,atoms%nat), overlap(0:atoms%ntype,atoms%ntype)
47 : REAL :: dispLength, maxDisp, limitDisp
48 : INTEGER :: iType, ierr, numDispReduce
49 : LOGICAL :: l_conv
50 :
51 : CHARACTER(len=100):: filename_add
52 :
53 : ! To calculate the current displacement
54 : TYPE(t_xml) :: xml
55 4 : TYPE(t_atoms) :: atoms_non_displaced
56 4 : TYPE(t_atoms) :: tempAtoms
57 :
58 4 : IF (fmpi%irank==0) THEN
59 2 : filename_add = ""
60 2 : IF (judft_was_argument("-add_name")) filename_add = TRIM(judft_string_for_argument("-add_name"))//"_"
61 2 : CALL xml%init(filename_add)
62 6 : ALLOCATE(pos(3,atoms%ntype,1));
63 :
64 6 : DO iType = 1, atoms%ntype
65 18 : pos(:,iType,1)=atoms%pos(:,atoms%firstAtom(iType))
66 : END DO
67 :
68 20 : ALLOCATE(force(3,atoms%ntype,1)); force(:,:,1)=force_new
69 2 : ALLOCATE(energies(1));energies(1)=energies_new
70 6 : ALLOCATE(displace(3,atoms%ntype),old_displace(3,atoms%ntype))
71 6 : ALLOCATE(tempDisplace(3,atoms%ntype),totalDisplace(3,atoms%ntype))
72 18 : totalDisplace=0.0
73 : ! Remove force components that are not selected for relaxation
74 6 : DO iType = 1, atoms%ntype
75 6 : IF (atoms%l_geo(iType)) THEN
76 16 : force(:,iType,1)=force(:,iType,1)*REAL(atoms%relax(:,iType))
77 : ELSE
78 0 : force(:,iType,1)=0.0
79 : END IF
80 : END DO
81 :
82 : ! Add history
83 2 : CALL read_relax(pos,force,energies)
84 :
85 : ! Determine new positions
86 2 : IF (SIZE(energies)==1.OR.input%forcemix==0) THEN
87 : ! No history present simple step
88 : ! choose a reasonable first guess for scaling
89 : ! this choice is based on a Debye temperature of 330K;
90 : ! modify as needed
91 : !alpha = (250.0/(MAXVAL(atoms%zatom)*input%xa))*((330./input%thetad)**2)
92 2 : CALL simple_step(input%forcealpha,0.25,force,displace)
93 0 : ELSE IF (input%forcemix==1) THEN
94 0 : CALL simple_cg(pos,force,displace)
95 0 : ELSE IF (input%forcemix==2) THEN
96 0 : CALL simple_bfgs(pos,force,displace)
97 : ELSE
98 0 : CALL juDFT_error('Unknown mixing scheme for forces', calledby='relaxation')
99 : END IF
100 :
101 : ! Check for convergence of forces/displacements
102 2 : maxDisp = 0.0
103 2 : l_conv = .TRUE.
104 6 : DO iType = 1, atoms%ntype
105 16 : IF (DOT_PRODUCT(force(:,iType,SIZE(force,3)),force(:,iType,SIZE(force,3)))>input%epsforce**2) l_conv=.FALSE.
106 16 : dispLength = SQRT(DOT_PRODUCT(displace(:,iType),displace(:,iType)))
107 4 : maxDisp = MAX(maxDisp,dispLength)
108 6 : IF (dispLength>input%epsdisp) l_conv=.FALSE.
109 : END DO
110 :
111 : ! Limit the maximal displacement in a single force relaxation step to limitDisp = 0.2 a_0.
112 2 : limitDisp = 0.2
113 2 : IF(maxDisp.GT.limitDisp) THEN
114 0 : displace(:,:) = limitDisp*displace(:,:) / maxDisp
115 : END IF
116 :
117 : ! New displacements relative to positions in inp.xml
118 : !CALL read_displacements(atoms,old_displace)
119 2 : CALL atoms_non_displaced%read_xml(xml)
120 2 : CALL xml%freeResources()
121 2 : CALL atoms_non_displaced%init(cell)
122 :
123 6 : DO iType = 1, atoms%ntype
124 : old_displace(:,iType) = atoms%pos(:,atoms%firstAtom(iType)) - &
125 18 : atoms_non_displaced%pos(:,atoms%firstAtom(iType))
126 : END DO
127 :
128 6 : tempAtoms = atoms_non_displaced
129 86 : tempDisplace = MATMUL(cell%bmat,totalDisplace)/tpi_const
130 :
131 2 : CALL rotate_to_all_sites(tempDisplace,atoms_non_displaced,cell,sym,dispAll)
132 22 : tempAtoms%taual(:,:)=atoms_non_displaced%taual(:,:)+dispAll(:,:)
133 222 : tempAtoms%pos=MATMUL(cell%amat,tempAtoms%taual)
134 :
135 6 : numDispReduce = 0
136 18 : overlap=1.0
137 22 : DO WHILE(ANY(overlap.GT.0.0))
138 18 : overlap = 0.0
139 20 : totalDisplace=displace+old_displace
140 :
141 10 : tempAtoms = atoms_non_displaced
142 122 : tempDisplace = MATMUL(cell%bmat,totalDisplace)/tpi_const
143 :
144 2 : CALL rotate_to_all_sites(tempDisplace,atoms_non_displaced,cell,sym,dispAll)
145 22 : tempAtoms%taual(:,:)=atoms_non_displaced%taual(:,:)+dispAll(:,:)
146 222 : tempAtoms%pos=MATMUL(cell%amat,tempAtoms%taual)
147 2 : CALL chkmt(tempAtoms,input,vacuum,cell ,.TRUE.,overlap=overlap)
148 :
149 20 : IF (ANY(overlap.GT.0.0)) THEN
150 0 : numDispReduce = numDispReduce + 1
151 0 : IF (numDispReduce.GE.3) THEN
152 0 : CALL juDFT_warn("Strong MT spheres crash in structural relaxation")
153 : END IF
154 0 : displace(:,:) = 0.5 * displace(:,:)
155 0 : WRITE(oUnit,*) 'Automatically reducing atom displacements because MT spheres crash into each other!'
156 0 : WRITE(*,*) 'Automatically reducing atom displacements because MT spheres crash into each other!'
157 : ! TODO: Why two calls?
158 : END IF
159 : END DO
160 :
161 : ! Write relax file
162 2 : CALL write_relax(pos,force,energies,totalDisplace)
163 :
164 : ! Structure in xsf-format
165 2 : OPEN (55,file="struct-relax.xsf",status='replace')
166 2 : CALL xsf_WRITE_atoms(55,tempAtoms,input%film,cell%amat)
167 2 : CLOSE (55)
168 : END IF
169 :
170 : #ifdef CPP_MPI
171 4 : CALL MPI_BCAST(l_conv,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
172 : #endif
173 :
174 4 : IF (l_conv) THEN
175 0 : CALL judft_end("Structural relaxation: Done",fmpi%irank)
176 : ELSE
177 4 : CALL mixing_history_reset(fmpi)
178 4 : CALL judft_end("Structural relaxation: new displacements generated",fmpi%irank)
179 : END IF
180 :
181 4 : END SUBROUTINE relaxation
182 :
183 2 : SUBROUTINE simple_step(alpha,maxdisp,force,displace)
184 :
185 : IMPLICIT NONE
186 :
187 : REAL, INTENT(IN) :: alpha,maxdisp
188 : REAL, INTENT(IN) :: force(:,:,:)
189 : REAL, INTENT(OUT) :: displace(:,:)
190 :
191 : REAL :: corr
192 :
193 18 : displace = alpha*force(:,:,SIZE(force,3))
194 18 : corr=maxdisp/maxval(abs(displace))
195 2 : IF (corr<1.0) displace = corr*alpha*force(:,:,size(force,3))
196 :
197 2 : END SUBROUTINE simple_step
198 :
199 0 : SUBROUTINE simple_bfgs(pos,force,shift)
200 : !--------------------------------------------------------------------------
201 : ! Simple BFGS method to calculate shift out of old positions and forces
202 : !--------------------------------------------------------------------------
203 : USE m_constants
204 :
205 : REAL,INTENT(IN) :: pos(:,:,:),force(:,:,:)
206 : REAL,INTENT(OUT) :: shift(:,:)
207 :
208 : INTEGER :: n, i, j, hist_length, n_force
209 : REAL, ALLOCATABLE :: h(:,:)
210 0 : REAL, ALLOCATABLE :: p(:), y(:), v(:)
211 : REAL :: py, yy, gamma
212 :
213 0 : n_force=3*SIZE(pos,2)
214 0 : ALLOCATE(h(n_force,n_force))
215 0 : ALLOCATE(p(n_force),y(n_force),v(n_force))
216 :
217 : ! Calculate approx. Hessian
218 :
219 : ! Initialize H
220 0 : h = 0.0
221 :
222 0 : DO n = 1,n_force
223 0 : h(n,n) = 1.0
224 : END DO
225 :
226 : ! Loop over all iterations (including current)
227 0 : hist_length=SIZE(pos,3)
228 :
229 0 : DO n = 2,hist_length
230 : ! Differences
231 0 : p(:) = RESHAPE(pos(:,:,n)-pos(:,:,n-1),(/SIZE(p)/))
232 0 : y(:) = RESHAPE(force(:,:,n-1)-force(:,:,n),(/SIZE(p)/))
233 :
234 : ! Get necessary inner products and H|y>
235 0 : py = DOT_PRODUCT(p,y)
236 0 : v = MATMUL(y,h)
237 0 : yy = DOT_PRODUCT(y,v)
238 :
239 : ! Check that update will leave h positive definite:
240 0 : IF (py <= 0.0) THEN
241 0 : WRITE (oUnit,*) ' bfgs: <p|y> < 0'
242 0 : WRITE (oUnit,*) ' check convergence of forces'
243 : ! Starting over with initial hessian
244 0 : h = 0.0
245 0 : DO j = 1,n_force
246 0 : h(j,j) = 1.0
247 : END DO
248 : CYCLE
249 : ELSE
250 : ! Update h
251 0 : IF (n == 2) THEN
252 0 : gamma = py/yy
253 : ELSE
254 : gamma = 1.0
255 : END IF
256 :
257 0 : DO j = 1,n_force
258 0 : DO i = 1,n_force
259 0 : h(i,j) = (h(i,j) - (v(i)*p(j)+p(i)*v(j))/py)*gamma + (1.+gamma*yy/py)*p(i)*p(j)/py
260 : END DO
261 : END DO
262 : END IF
263 : END DO
264 :
265 0 : y(:) = RESHAPE(force(:,:,hist_length),(/SIZE(p)/))
266 0 : shift = RESHAPE(MATMUL(y,h),SHAPE(shift))
267 0 : END SUBROUTINE simple_bfgs
268 :
269 0 : SUBROUTINE simple_cg(pos,force,shift)
270 :
271 : REAL, INTENT(IN) :: pos(:,:,:),force(:,:,:)
272 : REAL, INTENT(OUT) :: shift(:,:)
273 :
274 0 : REAL :: f1(3,SIZE(pos,2)),f2(3,SIZE(pos,2))
275 0 : REAL :: dist(3,SIZE(pos,2))
276 : REAL :: eps
277 : INTEGER :: n_old, i, j
278 :
279 0 : eps = 1.0e-9
280 :
281 0 : n_old = SIZE(pos,3)-1
282 :
283 0 : dist(:,:) = pos(:,:,n_old+1)-pos(:,:,n_old)
284 :
285 0 : DO i = 1, SIZE(pos,2)
286 0 : DO j = 1, 3
287 0 : IF(ABS(dist(j,i)).LT.eps) dist(j,i) = eps ! To avoid calculation of 0.0/0.0 below.
288 : END DO
289 : END DO
290 :
291 0 : f1 = (force(:,:,n_old+1)-force(:,:,n_old))/dist
292 0 : f2 = force(:,:,n_old+1)-f1*pos(:,:,n_old+1)
293 0 : shift = -1.*f2/f1-force(:,:,n_old+1)
294 0 : END SUBROUTINE simple_cg
295 :
296 0 : END MODULE m_relaxation
|
