Branch data Line data Source code
1 : : !!****m* ABINIT/m_abi_mixing
2 : : !! NAME
3 : : !! m_abi_mixing
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 2008-2026 ABINIT group (XG, DC, GMR)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : :
22 : : module m_abi_mixing
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : use m_linalg_interfaces
28 : : use m_xmpi
29 : :
30 : : use iso_c_binding, only : c_int16_t
31 : :
32 : : use m_time, only : timab
33 : : use m_io_tools, only : open_file
34 : :
35 : : implicit none
36 : :
37 : : private
38 : : !!***
39 : :
40 : : integer, parameter, public :: ABI_MIXING_NONE = 0
41 : : integer, parameter, public :: ABI_MIXING_EIG = 1
42 : : integer, parameter, public :: ABI_MIXING_SIMPLE = 2
43 : : integer, parameter, public :: ABI_MIXING_ANDERSON = 3
44 : : integer, parameter, public :: ABI_MIXING_ANDERSON_2 = 4
45 : : integer, parameter, public :: ABI_MIXING_CG_ENERGY = 5
46 : : integer, parameter, public :: ABI_MIXING_CG_ENERGY_2 = 6
47 : : integer, parameter, public :: ABI_MIXING_PULAY = 7
48 : :
49 : : integer, parameter, public :: ABI_MIXING_POTENTIAL = 0
50 : : integer, parameter, public :: ABI_MIXING_DENSITY = 1
51 : :
52 : : integer, parameter, public :: ABI_MIXING_REAL_SPACE = 1
53 : : integer, parameter, public :: ABI_MIXING_FOURRIER_SPACE = 2
54 : :
55 : : integer, parameter, public :: ABI_MIXING_PULAY_STORAGE_FULL = 0
56 : : integer, parameter, public :: ABI_MIXING_PULAY_STORAGE_DELTA = 1
57 : :
58 : :
59 : : type, public :: abi_mixing_object
60 : : integer :: iscf,useextfpmd,use_rcpaw
61 : : integer :: nfft, nspden, kind, space
62 : :
63 : : logical :: useprec
64 : : integer :: mffmem
65 : : character(len = fnlen) :: diskCache
66 : : integer :: n_index, n_fftgr, n_pulayit, n_pawmix,n_rcpawmix
67 : : integer :: pulayhist_storage
68 : :
69 : : integer, dimension(:), pointer :: i_rhor, i_vtrial, i_vresid, i_vrespc
70 : : real(dp), dimension(:,:,:), pointer :: f_fftgr, f_atm
71 : : real(sp), dimension(:,:,:), pointer :: f_fftgr_sp
72 : : real(sp), dimension(:,:), pointer :: f_fftgr_trial_sp
73 : : integer(c_int16_t), dimension(:,:,:), pointer :: f_fftgr_delta_i2
74 : : real(dp), dimension(:), pointer :: f_fftgr_delta_scale
75 : : real(dp), dimension(:,:), pointer :: f_paw
76 : :
77 : : real(dp),dimension(:),pointer :: f_extfpmd
78 : : real(dp),dimension(:,:), pointer :: f_rcpaw
79 : :
80 : : ! Private
81 : : integer :: n_atom
82 : : real(dp), pointer :: xred(:,:), dtn_pc(:,:)
83 : : end type abi_mixing_object
84 : :
85 : : public :: abi_mixing_new
86 : : public :: abi_mixing_deallocate
87 : :
88 : : public :: abi_mixing_use_disk_cache
89 : : public :: abi_mixing_use_moving_atoms
90 : : public :: abi_mixing_copy_current_step
91 : :
92 : : public :: abi_mixing_eval_allocate
93 : : public :: abi_mixing_eval
94 : : public :: abi_mixing_eval_deallocate
95 : : !!***
96 : :
97 : : contains
98 : : !!***
99 : :
100 : :
101 : : !!****f* m_abi_mixing/init_
102 : : !! NAME
103 : : !! init_
104 : : !!
105 : : !! FUNCTION
106 : : !! Initialize the object
107 : : !!
108 : : !! SOURCE
109 : :
110 : 8115 : subroutine init_(mix)
111 : :
112 : : !Arguments ------------------------------------
113 : : !scalars
114 : : type(abi_mixing_object), intent(out) :: mix
115 : : ! *************************************************************************
116 : :
117 : : ! Default values.
118 : 8115 : mix%iscf = ABI_MIXING_NONE
119 : 8115 : mix%mffmem = 1
120 : 8115 : mix%n_index = 0
121 : 8115 : mix%n_fftgr = 0
122 : 8115 : mix%n_pulayit = 7
123 : 8115 : mix%pulayhist_storage = ABI_MIXING_PULAY_STORAGE_FULL
124 : 8115 : mix%n_pawmix = 0
125 : 8115 : mix%n_atom = 0
126 : 8115 : mix%space = 0
127 : 8115 : mix%useprec = .true.
128 : 8115 : mix%useextfpmd = 0
129 : 8115 : mix%use_rcpaw = 0
130 : 8115 : mix%n_rcpawmix = 0
131 : :
132 : 8115 : call nullify_(mix)
133 : :
134 : 8115 : end subroutine init_
135 : : !!***
136 : :
137 : : !!****f* m_abi_mixing/nullify
138 : : !! NAME
139 : : !! nullify_
140 : : !!
141 : : !! FUNCTION
142 : : !! Nullify the pointers
143 : : !!
144 : : !! SOURCE
145 : :
146 : 16230 : subroutine nullify_(mix)
147 : :
148 : : !Arguments ------------------------------------
149 : : !scalars
150 : : type(abi_mixing_object), intent(inout) :: mix
151 : : ! *************************************************************************
152 : :
153 : : ! Nullify internal pointers.
154 : 16230 : nullify(mix%i_rhor)
155 : 16230 : nullify(mix%i_vtrial)
156 : 16230 : nullify(mix%i_vresid)
157 : 16230 : nullify(mix%i_vrespc)
158 : 16230 : nullify(mix%f_fftgr)
159 : 16230 : nullify(mix%f_fftgr_sp)
160 : 16230 : nullify(mix%f_fftgr_trial_sp)
161 : 16230 : nullify(mix%f_fftgr_delta_i2)
162 : 16230 : nullify(mix%f_fftgr_delta_scale)
163 : 16230 : nullify(mix%f_atm)
164 : 16230 : nullify(mix%f_paw)
165 : 16230 : nullify(mix%f_extfpmd)
166 : 16230 : nullify(mix%f_rcpaw)
167 : :
168 : : end subroutine nullify_
169 : : !!***
170 : :
171 : : !!****f* m_abi_mixing/abi_mixing_new
172 : : !! NAME
173 : : !! abi_mixing_new
174 : : !!
175 : : !! FUNCTION
176 : : !!
177 : : !! INPUTS
178 : : !!
179 : : !! OUTPUT
180 : : !!
181 : : !! NOTES
182 : : !!
183 : : !! SOURCE
184 : :
185 : 8115 : subroutine abi_mixing_new(mix, iscf, kind, space, nfft, nspden, &
186 : : & npawmix, errid, errmess, npulayit, useprec, pulayhist_storage)
187 : :
188 : : !Arguments ------------------------------------
189 : : !scalars
190 : : type(abi_mixing_object), intent(out) :: mix
191 : : integer, intent(in) :: iscf, kind, space, nfft, nspden, npawmix
192 : : integer, intent(out) :: errid
193 : : character(len = 500), intent(out) :: errmess
194 : : integer, intent(in), optional :: npulayit
195 : : integer, intent(in), optional :: pulayhist_storage
196 : : logical, intent(in), optional :: useprec
197 : :
198 : : !Local variables-------------------------------
199 : : !scalars
200 : : integer :: ii !, i_stat
201 : : character(len = *), parameter :: subname = "abi_mixing_new"
202 : : ! *************************************************************************
203 : :
204 : : ! Set default values.
205 : 8115 : call init_(mix)
206 : :
207 : : ! Argument checkings.
208 [ - + ]: 8115 : if (kind /= ABI_MIXING_POTENTIAL .and. kind /= ABI_MIXING_DENSITY) then
209 : 0 : errid = AB7_ERROR_MIXING_ARG
210 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
211 : 0 : & ' abi_mixing_set_arrays: ERROR -',ch10,&
212 : 0 : & ' Mixing must be done on density or potential only.'
213 : 0 : return
214 : : end if
215 [ - + ]: 8115 : if (space /= ABI_MIXING_REAL_SPACE .and. space /= ABI_MIXING_FOURRIER_SPACE) then
216 : 0 : errid = AB7_ERROR_MIXING_ARG
217 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
218 : 0 : & ' abi_mixing_set_arrays: ERROR -',ch10,&
219 : 0 : & ' Mixing must be done in real or Fourrier space only.'
220 : 0 : return
221 : : end if
222 : : if (iscf /= ABI_MIXING_EIG .and. iscf /= ABI_MIXING_SIMPLE .and. &
223 : : & iscf /= ABI_MIXING_ANDERSON .and. &
224 : : & iscf /= ABI_MIXING_ANDERSON_2 .and. &
225 : : & iscf /= ABI_MIXING_CG_ENERGY .and. &
226 : : & iscf /= ABI_MIXING_PULAY .and. &
227 [ + + + + : 8115 : & iscf /= ABI_MIXING_CG_ENERGY_2 .and. &
- + ]
228 : : & iscf /= ABI_MIXING_NONE) then
229 : 0 : errid = AB7_ERROR_MIXING_ARG
230 : 0 : write(errmess, "(A,I0,A)") "Unknown mixing scheme (", iscf, ")."
231 : 0 : return
232 : : end if
233 : 8115 : errid = AB7_NO_ERROR
234 : :
235 : : ! Mandatory arguments.
236 : 8115 : mix%iscf = iscf
237 : 8115 : mix%kind = kind
238 : 8115 : mix%space = space
239 : 8115 : mix%nfft = nfft
240 : 8115 : mix%nspden = nspden
241 : 8115 : mix%n_pawmix = npawmix
242 : :
243 : : ! Optional arguments.
244 [ - + ]: 8115 : if (present(useprec)) mix%useprec = useprec
245 [ + - ]: 8115 : if (present(pulayhist_storage)) mix%pulayhist_storage = pulayhist_storage
246 [ - + ]: 8115 : if (mix%pulayhist_storage /= ABI_MIXING_PULAY_STORAGE_FULL .and. &
247 : : & mix%pulayhist_storage /= ABI_MIXING_PULAY_STORAGE_DELTA) then
248 : 0 : errid = AB7_ERROR_MIXING_ARG
249 : 0 : write(errmess, "(A,I0,A)") "Unknown Pulay history storage mode (", mix%pulayhist_storage, ")."
250 : 0 : return
251 : : end if
252 : :
253 : : ! Set-up internal dimensions.
254 : : !These arrays are needed only in the self-consistent case
255 : : if (iscf == ABI_MIXING_NONE) then
256 : : ! For iscf==0, one additional vector is needed.
257 : : ! The index 1 is attributed to the new residual potential.
258 : 73 : mix%n_fftgr=1 ; mix%n_index=1
259 : : else if (iscf == ABI_MIXING_EIG) then
260 : : ! For iscf==1, five additional vectors are needed.
261 : : ! The index 1 is attributed to the old trial potential,
262 : : ! The new residual potential, and the new
263 : : ! preconditioned residual potential receive now a temporary index
264 : : ! The indices number 4 and 5 are attributed to work vectors.
265 : 3 : mix%n_fftgr=5 ; mix%n_index=1
266 : : else if(iscf == ABI_MIXING_SIMPLE) then
267 : : ! For iscf==2, three additional vectors are needed.
268 : : ! The index number 1 is attributed to the old trial vector
269 : : ! The new residual potential, and the new preconditioned
270 : : ! residual potential, receive now a temporary index.
271 : 31 : mix%n_fftgr=3 ; mix%n_index=1
272 [ - + ]: 31 : if (.not. mix%useprec) mix%n_fftgr = 2
273 : : else if(iscf == ABI_MIXING_ANDERSON) then
274 : : ! For iscf==3 , four additional vectors are needed.
275 : : ! The index number 1 is attributed to the old trial vector
276 : : ! The new residual potential, and the new and old preconditioned
277 : : ! residual potential, receive now a temporary index.
278 : 16 : mix%n_fftgr=4 ; mix%n_index=2
279 [ - + ]: 16 : if (.not. mix%useprec) mix%n_fftgr = 3
280 : : else if (iscf == ABI_MIXING_ANDERSON_2) then
281 : : ! For iscf==4 , six additional vectors are needed.
282 : : ! The indices number 1 and 2 are attributed to two old trial vectors
283 : : ! The new residual potential, and the new and two old preconditioned
284 : : ! residual potentials, receive now a temporary index.
285 : 2 : mix%n_fftgr=6 ; mix%n_index=3
286 [ - + ]: 2 : if (.not. mix%useprec) mix%n_fftgr = 5
287 : : else if(iscf == ABI_MIXING_CG_ENERGY .or. iscf == ABI_MIXING_CG_ENERGY_2) then
288 : : ! For iscf==5 or 6, ten additional vectors are needed
289 : : ! The index number 1 is attributed to the old trial vector
290 : : ! The index number 6 is attributed to the search vector
291 : : ! Other indices are attributed now. Altogether ten vectors
292 : 30 : mix%n_fftgr=10 ; mix%n_index=3
293 : : else if(iscf == ABI_MIXING_PULAY) then
294 : : ! For iscf==7, lot of additional vectors are needed
295 : : ! The index number 1 is attributed to the old trial vector
296 : : ! The index number 2 is attributed to the old residual
297 : : ! The indices number 2 and 3 are attributed to two old precond. residuals
298 : : ! Other indices are attributed now.
299 [ + - ]: 7960 : if (present(npulayit)) mix%n_pulayit = npulayit
300 : 7960 : mix%n_fftgr=2+2*mix%n_pulayit ; mix%n_index=1+mix%n_pulayit
301 [ - + ]: 7960 : if (.not. mix%useprec) mix%n_fftgr = 1+2*mix%n_pulayit
302 [ + + - + ]: 7960 : if (mix%pulayhist_storage == ABI_MIXING_PULAY_STORAGE_DELTA .and. .not. mix%useprec) then
303 : 0 : errid = AB7_ERROR_MIXING_ARG
304 : 0 : write(errmess, '(4a)' ) ch10,&
305 : 0 : & ' abi_mixing_new: ERROR -',ch10,&
306 : 0 : & ' Delta-encoded Pulay history currently requires preconditioned Pulay storage.'
307 : 0 : return
308 : : end if
309 : : end if ! iscf cases
310 : :
311 : : ! Allocate new arrays.
312 : : !allocate(mix%i_rhor(mix%n_index), stat = i_stat)
313 : : !call memocc_abi(i_stat, mix%i_rhor, 'mix%i_rhor', subname)
314 [ - + + - : 24345 : ABI_MALLOC(mix%i_rhor,(mix%n_index))
+ - + - -
+ - + ]
315 [ + + ]: 72066 : mix%i_rhor(:)=0
316 : : !allocate(mix%i_vtrial(mix%n_index), stat = i_stat)
317 : : !call memocc_abi(i_stat, mix%i_vtrial, 'mix%i_vtrial', subname)
318 [ - + + - : 16230 : ABI_MALLOC(mix%i_vtrial,(mix%n_index))
- + - + ]
319 [ + + ]: 72066 : mix%i_vtrial(:)=0
320 : : !allocate(mix%i_vresid(mix%n_index), stat = i_stat)
321 : : !call memocc_abi(i_stat, mix%i_vresid, 'mix%i_vresid', subname)
322 [ - + + - : 16230 : ABI_MALLOC(mix%i_vresid,(mix%n_index))
- + - + ]
323 [ + + ]: 72066 : mix%i_vresid(:)=0
324 : : !allocate(mix%i_vrespc(mix%n_index), stat = i_stat)
325 : : !call memocc_abi(i_stat, mix%i_vrespc, 'mix%i_vrespc', subname)
326 [ - + + - : 16230 : ABI_MALLOC(mix%i_vrespc,(mix%n_index))
- + - + ]
327 [ + + ]: 72066 : mix%i_vrespc(:)=0
328 : :
329 : : ! Setup initial values.
330 : : if (iscf == ABI_MIXING_NONE) then
331 : 73 : mix%i_vresid(1)=1
332 : : else if (iscf == ABI_MIXING_EIG) then
333 : 3 : mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2 ; mix%i_vrespc(1)=3
334 : : else if(iscf == ABI_MIXING_SIMPLE) then
335 : 31 : mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2 ; mix%i_vrespc(1)=3
336 [ - + ]: 31 : if (.not. mix%useprec) mix%i_vrespc(1)=2
337 : : else if(iscf == ABI_MIXING_ANDERSON) then
338 : 16 : mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2
339 [ + - ]: 16 : if (mix%useprec) then
340 : 16 : mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=4
341 : : else
342 : 0 : mix%i_vrespc(1)=2 ; mix%i_vrespc(2)=3
343 : : end if
344 : : else if (iscf == ABI_MIXING_ANDERSON_2) then
345 : 2 : mix%i_vtrial(1)=1 ; mix%i_vtrial(2)=2
346 : 2 : mix%i_vresid(1)=3
347 [ + - ]: 2 : if (mix%useprec) then
348 : 2 : mix%i_vrespc(1)=4 ; mix%i_vrespc(2)=5 ; mix%i_vrespc(3)=6
349 : : else
350 : 0 : mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=4 ; mix%i_vrespc(3)=5
351 : : end if
352 : : else if(iscf == ABI_MIXING_CG_ENERGY .or. iscf == ABI_MIXING_CG_ENERGY_2) then
353 : 30 : mix%n_fftgr=10 ; mix%n_index=3
354 : 30 : mix%i_vtrial(1)=1
355 : 30 : mix%i_vresid(1)=2 ; mix%i_vresid(2)=4 ; mix%i_vresid(3)=7
356 : 30 : mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=5 ; mix%i_vrespc(3)=8
357 : 30 : mix%i_rhor(2)=9 ; mix%i_rhor(3)=10
358 : : else if(iscf == ABI_MIXING_PULAY) then
359 [ + + ]: 63716 : do ii=1,mix%n_pulayit
360 : 63716 : mix%i_vtrial(ii)=2*ii-1 ; mix%i_vrespc(ii)=2*ii
361 : : end do
362 : 7960 : mix%i_vrespc(mix%n_pulayit+1)=2*mix%n_pulayit+1
363 : 7960 : mix%i_vresid(1)=2*mix%n_pulayit+2
364 [ - + ]: 7960 : if (.not. mix%useprec) mix%i_vresid(1)=2
365 : : end if ! iscf cases
366 : :
367 : 8115 : end subroutine abi_mixing_new
368 : : !!***
369 : :
370 : : !!****f* m_abi_mixing/abi_mixing_use_disk_cache
371 : : !! NAME
372 : : !! abi_mixing_use_disk_cache
373 : : !!
374 : : !! FUNCTION
375 : : !!
376 : : !! INPUTS
377 : : !!
378 : : !! OUTPUT
379 : : !!
380 : : !! NOTES
381 : : !! Obsolete?
382 : : !!
383 : : !! SOURCE
384 : :
385 : 0 : subroutine abi_mixing_use_disk_cache(mix, fnametmp_fft)
386 : :
387 : : !Arguments ------------------------------------
388 : : !scalars
389 : : type(abi_mixing_object), intent(inout) :: mix
390 : : character(len = *), intent(in) :: fnametmp_fft
391 : : ! *************************************************************************
392 : :
393 [ # # ]: 0 : if (len(trim(fnametmp_fft)) > 0) then
394 : 0 : mix%mffmem = 0
395 : 0 : write(mix%diskCache, "(A)") fnametmp_fft
396 : : else
397 : 0 : mix%mffmem = 1
398 : : end if
399 : :
400 : 0 : end subroutine abi_mixing_use_disk_cache
401 : : !!***
402 : :
403 : :
404 : : !!****f* m_abi_mixing/abi_mixing_use_moving_atoms
405 : : !! NAME
406 : : !! abi_mixing_use_moving_atoms
407 : : !!
408 : : !! FUNCTION
409 : : !!
410 : : !! INPUTS
411 : : !!
412 : : !! OUTPUT
413 : : !!
414 : : !! SIDE EFFECTS
415 : : !!
416 : : !! NOTES
417 : : !!
418 : : !! SOURCE
419 : :
420 : 110 : subroutine abi_mixing_use_moving_atoms(mix, natom, xred, dtn_pc)
421 : :
422 : : !Arguments ------------------------------------
423 : : !scalars
424 : : type(abi_mixing_object), intent(inout) :: mix
425 : : integer, intent(in) :: natom
426 : : real(dp), intent(in), target :: dtn_pc(3, natom)
427 : : real(dp), intent(in), target :: xred(3, natom)
428 : :
429 : : ! *************************************************************************
430 : :
431 : 110 : mix%n_atom = natom
432 : 110 : mix%dtn_pc => dtn_pc
433 : 110 : mix%xred => xred
434 : :
435 : 110 : end subroutine abi_mixing_use_moving_atoms
436 : : !!***
437 : :
438 : :
439 : : !!****f* m_abi_mixing/abi_mixing_copy_current_step
440 : : !! NAME
441 : : !! abi_mixing_copy_current_step
442 : : !!
443 : : !! FUNCTION
444 : : !!
445 : : !! INPUTS
446 : : !!
447 : : !! OUTPUT
448 : : !!
449 : : !! SIDE EFFECTS
450 : : !!
451 : : !! NOTES
452 : : !!
453 : : !! SOURCE
454 : 55884 : subroutine abi_mixing_copy_current_step(mix, arr_resid, errid, errmess, &
455 : 55994 : & arr_respc, arr_paw_resid, arr_paw_respc, arr_atm, &
456 : 0 : & arr_extfpmd_resid,arr_extfpmd_respc,arr_rcpaw_resid,arr_rcpaw_respc)
457 : :
458 : : !Arguments ------------------------------------
459 : : !scalars
460 : : type(abi_mixing_object), intent(inout) :: mix
461 : : real(dp), intent(in) :: arr_resid(mix%space * mix%nfft, mix%nspden)
462 : : integer, intent(out) :: errid
463 : : character(len = 500), intent(out) :: errmess
464 : : real(dp), intent(in), optional :: arr_respc(mix%space * mix%nfft, mix%nspden)
465 : : real(dp), intent(in), optional :: arr_paw_resid(mix%n_pawmix), arr_paw_respc(mix%n_pawmix)
466 : : real(dp), intent(in), optional :: arr_atm(3, mix%n_atom)
467 : : real(dp),intent(in), optional :: arr_extfpmd_resid,arr_extfpmd_respc
468 : : real(dp),intent(in), optional :: arr_rcpaw_resid(mix%n_rcpawmix*mix%use_rcpaw)
469 : : real(dp),intent(in), optional :: arr_rcpaw_respc(mix%n_rcpawmix*mix%use_rcpaw)
470 : : ! *************************************************************************
471 : :
472 : :
473 [ + - - + ]: 55884 : if (mix%n_fftgr>0 .and. (.not. associated(mix%f_fftgr))) then
474 : 0 : errid = AB7_ERROR_MIXING_ARG
475 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
476 : 0 : & ' abi_mixing_set_arr_current_step: ERROR (1) -',ch10,&
477 : 0 : & ' Working arrays not yet allocated.'
478 : 0 : return
479 : : end if
480 [ + + - + ]: 55884 : if (mix%n_pawmix>0 .and. (.not. associated(mix%f_paw))) then
481 : 0 : errid = AB7_ERROR_MIXING_ARG
482 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
483 : 0 : & ' abi_mixing_set_arr_current_step: ERROR (2) -',ch10,&
484 : 0 : & ' Working arrays not yet allocated.'
485 : 0 : return
486 : : end if
487 [ + + - + ]: 55884 : if (mix%n_atom>0 .and. (.not. associated(mix%f_atm))) then
488 : 0 : errid = AB7_ERROR_MIXING_ARG
489 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
490 : 0 : & ' abi_mixing_set_arr_current_step: ERROR (3) -',ch10,&
491 : 0 : & ' Working arrays not yet allocated.'
492 : 0 : return
493 : : end if
494 [ + + - + ]: 55884 : if (mix%n_rcpawmix>0 .and. (.not. associated(mix%f_rcpaw))) then
495 : 0 : errid = AB7_ERROR_MIXING_ARG
496 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
497 : 0 : & ' abi_mixing_set_arr_current_step: ERROR (3) -',ch10,&
498 : 0 : & ' Working arrays not yet allocated.'
499 : 0 : return
500 : : end if
501 : 55884 : errid = AB7_NO_ERROR
502 : :
503 [ + - ]: 55884 : if (mix%n_fftgr>0) then
504 [ + + + + ]: 55884 : if (mix%iscf == ABI_MIXING_PULAY .and. &
505 : : & mix%pulayhist_storage /= ABI_MIXING_PULAY_STORAGE_FULL) then
506 [ + - + + : 73746 : if (mix%i_vresid(1)>0) mix%f_fftgr(:,:,2) = arr_resid(:,:)
+ + ]
507 [ + - + - : 73746 : if (present(arr_respc).and.mix%i_vrespc(1)>0) mix%f_fftgr(:,:,1) = arr_respc(:,:)
+ + + + ]
508 : : else
509 [ + - + + : 1337747113 : if (mix%i_vresid(1)>0) mix%f_fftgr(:,:,mix%i_vresid(1)) = arr_resid(:,:)
+ + ]
510 [ + - + + : 1306181243 : if (present(arr_respc).and.mix%i_vrespc(1)>0) mix%f_fftgr(:,:,mix%i_vrespc(1)) = arr_respc(:,:)
+ + + + ]
511 : : end if
512 : : end if
513 [ + + ]: 55884 : if (mix%n_pawmix>0) then
514 [ - + - - : 14911 : if (present(arr_paw_resid).and.mix%i_vresid(1)>0) mix%f_paw(:, mix%i_vresid(1)) = arr_paw_resid(:)
- - ]
515 [ - + - - : 14911 : if (present(arr_paw_respc).and.mix%i_vrespc(1)>0) mix%f_paw(:, mix%i_vrespc(1)) = arr_paw_respc(:)
- - ]
516 : : end if
517 [ + + ]: 55884 : if (mix%n_atom>0) then
518 [ + + + - : 1456 : if (present(arr_atm).and.mix%i_vresid(1)>0) mix%f_atm(:,:, mix%i_vresid(1)) = arr_atm(:,:)
+ + + + ]
519 : : end if
520 [ + - ]: 55884 : if(associated(mix%f_extfpmd)) then
521 [ - + - - ]: 55884 : if(present(arr_extfpmd_resid).and.mix%i_vresid(1)>0) mix%f_extfpmd(mix%i_vresid(1))=arr_extfpmd_resid
522 [ - + - - ]: 55884 : if(present(arr_extfpmd_respc).and.mix%i_vrespc(1)>0) mix%f_extfpmd(mix%i_vrespc(1))=arr_extfpmd_respc
523 : : endif
524 [ + + ]: 55884 : if(mix%n_rcpawmix>0) then
525 [ - + - - : 18 : if(present(arr_rcpaw_resid).and.mix%i_vresid(1)>0) mix%f_rcpaw(:,mix%i_vresid(1))=arr_rcpaw_resid(:)
- - ]
526 [ # # # # ]: 0 : if(present(arr_rcpaw_resid).and.mix%i_vrespc(1)>0)mix%f_rcpaw(:,mix%i_vrespc(1))=arr_rcpaw_respc(:)
527 : : endif
528 : :
529 [ + - - + : 111878 : end subroutine abi_mixing_copy_current_step
- + + + -
+ - + ]
530 : : !!***
531 : :
532 : :
533 : : !!****f* m_abi_mixing/abi_mixing_eval_allocate
534 : : !! NAME
535 : : !! abi_mixing_eval_allocate
536 : : !!
537 : : !! FUNCTION
538 : : !!
539 : : !! INPUTS
540 : : !!
541 : : !! OUTPUT
542 : : !!
543 : : !! SIDE EFFECTS
544 : : !!
545 : : !! NOTES
546 : : !!
547 : : !! SOURCE
548 : :
549 : 55884 : subroutine abi_mixing_eval_allocate(mix, istep)
550 : :
551 : : !Arguments ------------------------------------
552 : : !scalars
553 : : type(abi_mixing_object), intent(inout) :: mix
554 : : integer, intent(in), optional :: istep
555 : :
556 : : !Local variables-------------------------------
557 : : !scalars
558 : : integer :: istep_,temp_unit !, i_stat
559 : : real(dp) :: tsec(2)
560 : : character(len = *), parameter :: subname = "abi_mixing_eval_allocate"
561 : : character(len=500) :: msg
562 : :
563 : : ! *************************************************************************
564 : :
565 : 55884 : istep_ = 1
566 [ + - ]: 55884 : if (present(istep)) istep_ = istep
567 : :
568 : : ! Allocate work array.
569 [ + + ]: 55884 : if (.not. associated(mix%f_fftgr)) then
570 : : !allocate(mix%f_fftgr(mix%space * mix%nfft,mix%nspden,mix%n_fftgr), stat = i_stat)
571 : : !call memocc_abi(i_stat, mix%f_fftgr, 'mix%f_fftgr', subname)
572 [ + + + + ]: 7853 : if (mix%iscf == ABI_MIXING_PULAY .and. &
573 : : & mix%pulayhist_storage /= ABI_MIXING_PULAY_STORAGE_FULL) then
574 [ - + + - : 5 : ABI_MALLOC(mix%f_fftgr,(mix%space * mix%nfft,mix%nspden,2))
- + + - +
- + - + -
- + - + ]
575 : : else
576 [ - + + - : 39260 : ABI_MALLOC(mix%f_fftgr,(mix%space * mix%nfft,mix%nspden,mix%n_fftgr))
- + + - -
+ + - + -
+ - - + -
+ ]
577 : : end if
578 [ + + + + : 2794335936 : mix%f_fftgr(:,:,:)=zero
+ + ]
579 [ - + - - : 7853 : if (mix%mffmem == 0 .and. istep_ > 1 .and. mix%n_fftgr>0) then
- - ]
580 : 0 : call timab(83,1,tsec)
581 [ # # ]: 0 : if (open_file(mix%diskCache,msg,newunit=temp_unit,form='unformatted',status='old') /= 0) then
582 : 0 : ABI_ERROR(msg)
583 : : end if
584 : 0 : rewind(temp_unit)
585 [ # # # # : 0 : read(temp_unit) mix%f_fftgr
# # ]
586 [ # # ]: 0 : if (mix%n_pawmix == 0) close(unit=temp_unit)
587 : 0 : call timab(83,2,tsec)
588 : : end if
589 : : end if
590 [ + + + + ]: 55884 : if (mix%iscf == ABI_MIXING_PULAY .and. &
591 : : & mix%pulayhist_storage == ABI_MIXING_PULAY_STORAGE_DELTA) then
592 [ + + ]: 9 : if (.not. associated(mix%f_fftgr_sp)) then
593 [ - + + - : 5 : ABI_MALLOC(mix%f_fftgr_sp,(mix%space * mix%nfft,mix%nspden,mix%n_pulayit+1))
- + + - +
- - + + -
+ - + - -
+ - + ]
594 [ + + + + : 65553 : mix%f_fftgr_sp(:,:,:)=zero_sp
+ + ]
595 : : end if
596 [ + + ]: 9 : if (.not. associated(mix%f_fftgr_trial_sp)) then
597 [ - + + - : 4 : ABI_MALLOC(mix%f_fftgr_trial_sp,(mix%space * mix%nfft,mix%nspden))
- + + - +
- + - - +
- + ]
598 [ + + + + ]: 8194 : mix%f_fftgr_trial_sp(:,:)=zero_sp
599 : : end if
600 [ + + ]: 9 : if (.not. associated(mix%f_fftgr_delta_i2)) then
601 [ - + + - : 5 : ABI_MALLOC(mix%f_fftgr_delta_i2,(mix%space * mix%nfft,mix%nspden,mix%n_pulayit))
- + + - -
+ + - + -
+ - - + -
+ ]
602 [ + + + + : 57359 : mix%f_fftgr_delta_i2(:,:,:)=0_c_int16_t
+ + ]
603 : : end if
604 [ + + ]: 9 : if (.not. associated(mix%f_fftgr_delta_scale)) then
605 [ - + + - : 3 : ABI_MALLOC(mix%f_fftgr_delta_scale,(mix%n_pulayit))
+ - + - -
+ - + ]
606 [ + + ]: 8 : mix%f_fftgr_delta_scale(:)=one
607 : : end if
608 : : end if
609 : : ! Allocate PAW work array.
610 [ + + ]: 55884 : if (.not. associated(mix%f_paw)) then
611 : : !allocate(mix%f_paw(mix%n_pawmix,mix%n_fftgr), stat = i_stat)
612 : : !call memocc_abi(i_stat, mix%f_paw, 'mix%f_paw', subname)
613 [ + + + - : 31412 : ABI_MALLOC(mix%f_paw,(mix%n_pawmix,mix%n_fftgr))
- + + - +
- + + - +
- + ]
614 [ + + + - ]: 7853 : if (mix%n_pawmix > 0 .and. mix%n_fftgr>0) then
615 [ + + + + ]: 5439361 : mix%f_paw(:,:)=zero
616 [ - + - - ]: 1517 : if (mix%mffmem == 0 .and. istep_ > 1) then
617 [ # # # # ]: 0 : read(temp_unit) mix%f_paw
618 : 0 : close(unit=temp_unit)
619 : 0 : call timab(83,2,tsec)
620 : : end if
621 : : end if
622 : : end if
623 : : ! Allocate atom work array.
624 [ + + ]: 55884 : if (.not. associated(mix%f_atm)) then
625 : : !allocate(mix%f_atm(3,mix%n_atom,mix%n_fftgr), stat = i_stat)
626 : : !call memocc_abi(i_stat, mix%f_atm, 'mix%f_atm', subname)
627 [ + + + - : 31412 : ABI_MALLOC(mix%f_atm,(3,mix%n_atom,mix%n_fftgr))
- + + - +
- + + - +
- + ]
628 : : end if
629 : :
630 [ + + ]: 55884 : if(.not.associated(mix%f_extfpmd)) then
631 [ + + ]: 7853 : if(mix%useextfpmd==1) then
632 [ - + + - : 21 : ABI_MALLOC(mix%f_extfpmd,(mix%n_fftgr))
+ - + - -
+ - + ]
633 [ + + ]: 119 : mix%f_extfpmd=zero
634 : : else
635 [ - + ]: 7846 : ABI_MALLOC(mix%f_extfpmd,(0))
636 : : endif
637 : : endif
638 : :
639 [ + + ]: 55884 : if(.not.associated(mix%f_rcpaw)) then
640 [ + + ]: 7853 : if(mix%use_rcpaw==1) then
641 [ - + + - : 12 : ABI_MALLOC(mix%f_rcpaw,(mix%n_rcpawmix,mix%n_fftgr))
- + + - +
- + - - +
- + ]
642 [ + + + + ]: 195 : mix%f_rcpaw=zero
643 : : else
644 [ - + ]: 7850 : ABI_MALLOC(mix%f_rcpaw,(0,0))
645 : : endif
646 : : endif
647 : :
648 : 55884 : end subroutine abi_mixing_eval_allocate
649 : : !!***
650 : :
651 : :
652 : : !!****f* m_abi_mixing/abi_mixing_eval_deallocate
653 : : !! NAME
654 : : !! abi_mixing_eval_deallocate
655 : : !!
656 : : !! FUNCTION
657 : : !!
658 : : !! INPUTS
659 : : !!
660 : : !! OUTPUT
661 : : !!
662 : : !! SIDE EFFECTS
663 : : !!
664 : : !! NOTES
665 : : !!
666 : : !! SOURCE
667 : :
668 : 55884 : subroutine abi_mixing_eval_deallocate(mix)
669 : :
670 : : !Arguments ------------------------------------
671 : : !scalars
672 : : type(abi_mixing_object), intent(inout) :: mix
673 : :
674 : : !Local variables-------------------------------
675 : : !scalars
676 : : integer :: temp_unit !i_all, i_stat
677 : : real(dp) :: tsec(2)
678 : : character(len = *), parameter :: subname = "abi_mixing_eval_deallocate"
679 : : character(len=500) :: msg
680 : :
681 : : ! *************************************************************************
682 : :
683 : : ! Save on disk and deallocate work array in case on disk cache only.
684 [ - + ]: 55884 : if (mix%mffmem == 0) then
685 : 0 : call timab(83,1,tsec)
686 [ # # ]: 0 : if (open_file(mix%diskCache,msg,newunit=temp_unit,form='unformatted',status='unknown') /= 0) then
687 : 0 : ABI_ERROR(msg)
688 : : end if
689 : 0 : rewind(temp_unit)
690 : : ! VALGRIND complains not all of f_fftgr_disk is initialized
691 [ # # ]: 0 : if (mix%n_fftgr > 0) then
692 [ # # # # : 0 : write(temp_unit) mix%f_fftgr
# # ]
693 : : end if
694 [ # # # # ]: 0 : if (mix%n_pawmix > 0 .and. mix%n_fftgr > 0) then
695 [ # # # # ]: 0 : write(temp_unit) mix%f_paw
696 : : end if
697 : 0 : close(unit=temp_unit)
698 : 0 : call timab(83,2,tsec)
699 [ # # ]: 0 : if (associated(mix%f_fftgr)) then
700 : 0 : ABI_FREE(mix%f_fftgr)
701 : : nullify(mix%f_fftgr)
702 : : end if
703 [ # # ]: 0 : if (associated(mix%f_paw)) then
704 : 0 : ABI_FREE(mix%f_paw)
705 : : nullify(mix%f_paw)
706 : : end if
707 [ # # ]: 0 : if(associated(mix%f_extfpmd)) then
708 : 0 : ABI_FREE(mix%f_extfpmd)
709 : : nullify(mix%f_extfpmd)
710 : : endif
711 [ # # ]: 0 : if(associated(mix%f_rcpaw)) then
712 : 0 : ABI_FREE(mix%f_rcpaw)
713 : : nullify(mix%f_rcpaw)
714 : : endif
715 : : end if
716 : :
717 : 55884 : end subroutine abi_mixing_eval_deallocate
718 : : !!***
719 : :
720 : :
721 : : !!****f* m_abi_mixing/abi_mixing_eval
722 : : !! NAME
723 : : !! abi_mixing_eval
724 : : !!
725 : : !! FUNCTION
726 : : !!
727 : : !! INPUTS
728 : : !!
729 : : !! OUTPUT
730 : : !!
731 : : !! SIDE EFFECTS
732 : : !!
733 : : !! NOTES
734 : : !!
735 : : !! SOURCE
736 : :
737 : 111768 : subroutine abi_mixing_eval(mix, arr, istep, nfftot, ucvol, &
738 : : & mpi_comm, mpi_summarize, errid, errmess, &
739 : 55281 : & reset, isecur, pawarr, pawopt, response, etotal, potden, &
740 : 13329 : & resnrm, nelect_extfpmd,rcpaw_arr,comm_atom)
741 : :
742 : : !Arguments ------------------------------------
743 : : !scalars
744 : : type(abi_mixing_object), intent(inout) :: mix
745 : : integer, intent(in) :: istep, nfftot, mpi_comm
746 : : real(dp), intent(in) :: ucvol
747 : : real(dp), intent(inout) :: arr(mix%space * mix%nfft,mix%nspden)
748 : : logical, intent(in) :: mpi_summarize
749 : : integer, intent(out) :: errid
750 : : character(len = 500), intent(out) :: errmess
751 : : logical, intent(in), optional :: reset
752 : : integer, intent(in), optional :: isecur, comm_atom, pawopt, response
753 : : real(dp), intent(inout), optional, target :: pawarr(mix%n_pawmix)
754 : : real(dp), intent(in), optional :: etotal
755 : : real(dp), intent(in), optional :: potden(mix%space * mix%nfft,mix%nspden)
756 : : real(dp), intent(out), optional :: resnrm
757 : : real(dp),intent(inout),optional :: nelect_extfpmd
758 : : real(dp),intent(inout),optional,target :: rcpaw_arr(mix%n_rcpawmix*mix%use_rcpaw)
759 : :
760 : : !Local variables-------------------------------
761 : : !scalars
762 : : integer :: moveAtm, dbl_nnsclo, initialized, isecur_, comm_atom_
763 : : integer :: usepaw, pawoptmix_, response_
764 : : real(dp) :: resnrm_,nelect_extfpmd_
765 : : !arrays
766 : : real(dp),target :: dum(1),dum0(0)
767 : 55884 : real(dp),pointer :: pawarr_(:),rcpawarr_(:)
768 : :
769 : : ! *************************************************************************
770 : :
771 : : ! Argument checkings.
772 : : !if (mix%iscf == ABI_MIXING_NONE) then
773 : : ! errid = AB7_ERROR_MIXING_ARG
774 : : ! write(errmess, '(a,a,a,a)' )ch10,&
775 : : ! & ' abi_mixing_eval: ERROR -',ch10,&
776 : : ! & ' No method has been chosen.'
777 : : ! return
778 : : !end if
779 [ - + ]: 14911 : if (mix%n_pawmix > 0 .and. .not. present(pawarr)) then
780 : 0 : errid = AB7_ERROR_MIXING_ARG
781 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
782 : 0 : & ' abi_mixing_eval: ERROR -',ch10,&
783 : 0 : & ' PAW is used, but no pawarr argument provided.'
784 : 0 : return
785 : : end if
786 [ + + + - : 55884 : if (mix%n_atom > 0 .and. (.not. associated(mix%dtn_pc) .or. .not. associated(mix%xred))) then
- + ]
787 : 0 : errid = AB7_ERROR_MIXING_ARG
788 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
789 : 0 : & ' abi_mixing_eval: ERROR -',ch10,&
790 : 0 : & ' Moving atoms is used, but no xred or dtn_pc attributes provided.'
791 : 0 : return
792 : : end if
793 : 55884 : errid = AB7_NO_ERROR
794 : :
795 : : ! Reset if requested
796 : 55884 : initialized = 1
797 [ + - ]: 55884 : if (present(reset)) then
798 [ + + ]: 55884 : if (reset) initialized = 0
799 : : end if
800 : :
801 : : ! Miscellaneous
802 : 55884 : moveAtm = 0
803 [ + + ]: 55884 : if (mix%n_atom > 0) moveAtm = 1
804 : 55884 : isecur_ = 0
805 [ + + ]: 55884 : if (present(isecur)) isecur_ = isecur
806 : 55884 : comm_atom_ = xmpi_comm_self
807 [ + + ]: 55884 : if (present(comm_atom)) comm_atom_ = comm_atom
808 : 55884 : usepaw = 0
809 [ + + ]: 55884 : if (mix%n_pawmix > 0) usepaw = 1
810 : 55884 : pawoptmix_ = 0
811 [ + + ]: 55884 : if (present(pawopt)) pawoptmix_ = pawopt
812 : 55884 : response_ = 0
813 [ + + ]: 55884 : if (present(response)) response_ = response
814 [ + + ]: 55884 : pawarr_ => dum ; if (present(pawarr)) pawarr_ => pawarr
815 [ + + ]: 55884 : nelect_extfpmd_=zero;if(present(nelect_extfpmd)) nelect_extfpmd_=nelect_extfpmd
816 [ + + ]: 55884 : rcpawarr_ => dum0; if(present(rcpaw_arr)) rcpawarr_ => rcpaw_arr
817 : :
818 : : ! Do the mixing.
819 : 55884 : resnrm_ = 0.d0
820 [ + + ]: 55884 : if (mix%iscf == ABI_MIXING_NONE) then
821 [ + + + + ]: 31566473 : arr(:,:)=arr(:,:)+mix%f_fftgr(:,:,1)
822 [ + + ]: 55281 : else if (mix%iscf == ABI_MIXING_EIG) then
823 : : ! This routine compute the eigenvalues of the SCF operator
824 : : call scfeig(istep, mix%space * mix%nfft, mix%nspden, &
825 : : & mix%f_fftgr(:,:,mix%i_vrespc(1)), arr, &
826 [ - + - - : 47 : & mix%f_fftgr(:,:,1), mix%f_fftgr(:,:,4:5), errid, errmess)
- - - - -
+ - - - -
- - - + -
- - - - -
- - - + -
- - - - +
- - - - -
+ - - - -
- - ]
827 : : else if (mix%iscf == ABI_MIXING_SIMPLE .or. &
828 : : & mix%iscf == ABI_MIXING_ANDERSON .or. &
829 [ + + ]: 55234 : & mix%iscf == ABI_MIXING_ANDERSON_2 .or. &
830 : : & mix%iscf == ABI_MIXING_PULAY) then
831 [ + + + + ]: 55048 : if (mix%iscf == ABI_MIXING_PULAY .and. &
832 : : & mix%pulayhist_storage == ABI_MIXING_PULAY_STORAGE_DELTA) then
833 : : call scfopt_pulay_delta_sp(mix%space, mix%f_fftgr, mix%f_fftgr_sp, &
834 : : & mix%f_fftgr_trial_sp, mix%f_fftgr_delta_i2, mix%f_fftgr_delta_scale, &
835 : : & mix%f_paw, istep, mix%i_vrespc, mix%i_vtrial, mpi_comm, &
836 : : & mpi_summarize, mix%nfft, mix%n_pawmix, mix%nspden, mix%n_fftgr, &
837 : : & mix%n_index, mix%kind, pawoptmix_, usepaw, pawarr_, resnrm_, &
838 : : & arr, errid, errmess, mix%useextfpmd, mix%f_extfpmd, &
839 : : & nelect_extfpmd_, mix%use_rcpaw, mix%n_rcpawmix, mix%f_rcpaw, &
840 [ + - - - : 9 : & rcpawarr_, comm_atom_)
- - - - +
- - - - -
+ - - - -
- + - - -
- - + - -
- - - - -
- + - - -
- - + - -
- + - - -
+ - + - -
- + - - -
- - + ]
841 [ + - ]: 55039 : else if (present(comm_atom)) then
842 : : call scfopt(mix%space, mix%f_fftgr,mix%f_paw,mix%iscf,istep,&
843 : : & mix%i_vrespc,mix%i_vtrial, &
844 : : & mpi_comm,mpi_summarize,mix%nfft,mix%n_pawmix,mix%nspden, &
845 : : & mix%n_fftgr,mix%n_index,mix%kind,pawoptmix_,usepaw,pawarr_, &
846 : : & resnrm_, arr, errid, errmess,&
847 : : & mix%useextfpmd,mix%f_extfpmd,nelect_extfpmd_,&
848 : : & mix%use_rcpaw,mix%n_rcpawmix,mix%f_rcpaw,rcpawarr_,&
849 [ + - - - : 55039 : & comm_atom=comm_atom)
- - - - -
- + - - -
- - - - +
- - - - -
+ - - - -
- + - - -
- - + - -
- - - - -
- + - - -
- - - - +
- - - - -
+ - - - +
- - - + -
+ - - - +
- - - - -
+ ]
850 : : else
851 : : call scfopt(mix%space, mix%f_fftgr,mix%f_paw,mix%iscf,istep,&
852 : : & mix%i_vrespc,mix%i_vtrial, &
853 : : & mpi_comm,mpi_summarize,mix%nfft,mix%n_pawmix,mix%nspden, &
854 : : & mix%n_fftgr,mix%n_index,mix%kind,pawoptmix_,usepaw,pawarr_, &
855 : : & resnrm_, arr, errid, errmess,&
856 : : & mix%useextfpmd,mix%f_extfpmd,nelect_extfpmd_,&
857 [ # # # # : 0 : & mix%use_rcpaw,mix%n_rcpawmix,mix%f_rcpaw,rcpawarr_)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
858 : : end if
859 : : ! Change atomic positions
860 [ + + + + : 55048 : if((istep==1 .or. mix%iscf==ABI_MIXING_SIMPLE) .and. mix%n_atom > 0)then
+ + ]
861 : : ! GAF: 2009-06-03
862 : : ! Apparently there are not reason
863 : : ! to restrict iscf=2 for ionmov=5
864 [ + + + + ]: 1266 : mix%xred(:,:) = mix%xred(:,:) + mix%dtn_pc(:,:)
865 : : end if
866 [ + - ]: 186 : else if (mix%iscf == ABI_MIXING_CG_ENERGY .or. mix%iscf == ABI_MIXING_CG_ENERGY_2) then
867 : : ! Optimize next vtrial using an algorithm based
868 : : ! on the conjugate gradient minimization of etotal
869 [ - + ]: 186 : if (.not. present(etotal) .or. .not. present(potden)) then
870 : 0 : errid = AB7_ERROR_MIXING_ARG
871 : 0 : write(errmess, '(a,a,a,a)' )ch10,&
872 : 0 : & ' abi_mixing_eval: ERROR -',ch10,&
873 : 0 : & ' Arguments etotal or potden are missing for CG on energy methods.'
874 : 0 : return
875 : : end if
876 [ + + ]: 186 : if (mix%n_atom == 0) then
877 [ - + ]: 166 : ABI_MALLOC(mix%xred,(3,0))
878 [ - + ]: 166 : ABI_MALLOC(mix%dtn_pc,(3,0))
879 : : end if
880 : : call scfcge(mix%space,dbl_nnsclo,mix%dtn_pc,etotal,mix%f_atm,&
881 : : & mix%f_fftgr,initialized,mix%iscf,isecur_,istep,&
882 : : & mix%i_rhor,mix%i_vresid,mix%i_vrespc,moveAtm,&
883 : : & mpi_comm,mpi_summarize,mix%n_atom,mix%nfft,nfftot,&
884 : : & mix%nspden,mix%n_fftgr,mix%n_index,mix%kind,&
885 [ + - - - : 186 : & response_,potden,ucvol,arr,mix%xred, errid, errmess)
- - - - +
- - - - -
- - - - +
- - - - -
- - - - +
- - - - -
+ - - - -
- + - - -
- - + - -
- - - - -
- + - + -
- - - - -
- + - - -
- - - - +
- - - + -
- - + - -
- + - - -
- ]
886 [ + + ]: 186 : if (mix%n_atom == 0) then
887 [ - + ]: 166 : ABI_FREE(mix%xred)
888 [ - + ]: 166 : ABI_FREE(mix%dtn_pc)
889 : : end if
890 [ - + ]: 186 : if (dbl_nnsclo == 1) errid = AB7_ERROR_MIXING_INC_NNSLOOP
891 : : end if
892 : :
893 [ + + ]: 55884 : if(present(nelect_extfpmd)) nelect_extfpmd=nelect_extfpmd_
894 [ - + ]: 55884 : if (present(resnrm)) resnrm = resnrm_
895 : : nullify(rcpawarr_)
896 : :
897 [ + + + + : 124494 : end subroutine abi_mixing_eval
+ + ]
898 : : !!***
899 : :
900 : :
901 : : !!****f* m_abi_mixing/abi_mixing_deallocate
902 : : !! NAME
903 : : !! abi_mixing_deallocate
904 : : !!
905 : : !! FUNCTION
906 : : !!
907 : : !! INPUTS
908 : : !!
909 : : !! OUTPUT
910 : : !!
911 : : !! SIDE EFFECTS
912 : : !!
913 : : !! NOTES
914 : : !!
915 : : !! SOURCE
916 : :
917 : 8115 : subroutine abi_mixing_deallocate(mix)
918 : :
919 : : !Arguments ------------------------------------
920 : : !scalars
921 : : type(abi_mixing_object), intent(inout) :: mix
922 : :
923 : : !Local variables-------------------------------
924 : : !scalars
925 : : character(len = *), parameter :: subname = "abi_mixing_deallocate"
926 : : ! *************************************************************************
927 : :
928 [ + - ]: 8115 : ABI_SFREE_PTR(mix%i_rhor)
929 [ + - ]: 8115 : ABI_SFREE_PTR(mix%i_vtrial)
930 [ + - ]: 8115 : ABI_SFREE_PTR(mix%i_vresid)
931 [ + - ]: 8115 : ABI_SFREE_PTR(mix%i_vrespc)
932 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_fftgr)
933 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_fftgr_sp)
934 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_fftgr_trial_sp)
935 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_fftgr_delta_i2)
936 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_fftgr_delta_scale)
937 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_paw)
938 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_atm)
939 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_extfpmd)
940 [ + + ]: 8115 : ABI_SFREE_PTR(mix%f_rcpaw)
941 : :
942 : 8115 : call nullify_(mix)
943 : :
944 : 8115 : end subroutine abi_mixing_deallocate
945 : : !!***
946 : :
947 : : !!****f* m_abi_mixing/scfcge
948 : : !!
949 : : !! NAME
950 : : !! scfcge
951 : : !!
952 : : !! FUNCTION
953 : : !! Compute the next vtrial of the SCF cycle.
954 : : !! Uses a conjugate gradient minimization of the total energy
955 : : !! Can move only the trial potential (if moved_atm_inside==0), or
956 : : !! move the trial atomic positions as well (if moved_atm_inside==1).
957 : : !!
958 : : !! INPUTS
959 : : !! cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
960 : : !! dtn_pc(3,natom)=preconditioned change of atomic position, in reduced
961 : : !! coordinates. Will be quickly transferred to f_atm(:,:,i_vrespc(1))
962 : : !! etotal=the actual total energy
963 : : !! initialized= if 0, the initialization of the gstate run is not yet finished
964 : : !! iscf =5 => SCF cycle, CG based on estimation of energy gradient
965 : : !! =6 => SCF cycle, CG based on true minimization of the energy
966 : : !! isecur=level of security of the computation
967 : : !! istep= number of the step in the SCF cycle
968 : : !! moved_atm_inside: if==1, the atoms are allowed to move.
969 : : !! mpicomm=the mpi communicator used for the summation
970 : : !! mpi_summarize=set it to .true. if parallelisation is done over FFT
971 : : !! natom=number of atoms
972 : : !! nfft=(effective) number of FFT grid points (for this processor)
973 : : !! nfftot=total number of FFT grid points
974 : : !! nspden=number of spin-density components
975 : : !! n_fftgr=third dimension of the array f_fftgr
976 : : !! n_index=dimension for indices of potential/density (see i_vresid, ivrespc, i_rhor...)
977 : : !! opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential
978 : : !! 1 vtrial (and also f_fftgr) actually contains the trial density
979 : : !! response= if 0, GS calculation, if 1, RF calculation, intrinsically harmonic !
980 : : !! rhor(cplex*nfft,nspden)=actual density
981 : : !! ucvol=unit cell volume in bohr**3
982 : : !!
983 : : !! OUTPUT
984 : : !! dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.
985 : : !!
986 : : !! SIDE EFFECTS
987 : : !! Input/Output:
988 : : !! vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave
989 : : !! the input residual of the potential and Hellman-Feynman forces
990 : : !! at output, it is the new trial potential .
991 : : !! xred(3,natom)=(needed if moved_atm_inside==1)
992 : : !! reduced dimensionless atomic coordinates
993 : : !! at input, those that generated the input residual of the potential
994 : : !! and Hellman-Feynman forces, at output, these are the new ones.
995 : : !! f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid :
996 : : !! The input vtrial is transferred, at output, in f_fftgr(:,:,1).
997 : : !! The input f_fftgr(:,:,i_vresid(1)) contains the last residual.
998 : : !! the value of i_vresid(1) is transferred to i_vresid(2) at output.
999 : : !! The input f_fftgr(:,:,i_vresid(2)) contains the old residual.
1000 : : !! the value of i_vresid(2) is transferred to i_vresid(3) at output.
1001 : : !! The input f_fftgr(:,:,i_vresid(3)) contains the previous last residual.
1002 : : !! For the preconditioned potential residual, the same logic as for the
1003 : : !! the potential residual is used, with i_vrespc replacing i_vresid.
1004 : : !! The input rhor is transferred, at output, in f_fft(:,:,i_rhor(2)).
1005 : : !! The old density is input in f_fft(:,:,i_rhor(2)), and the value of
1006 : : !! i_rhor(2) is transferred to i_rhor(3) before the end of the routine.
1007 : : !! The input/output search vector is stored in f_fftgr(:,:,6)
1008 : : !! f_atm(3,natom,n_fftgr)=different functions defined for each atom :
1009 : : !! The input xred is transferred, at output, in f_atm(:,:,1).
1010 : : !! The input f_atm(:,:,i_vresid(1)) contains minus the HF forces.
1011 : : !! the value of i_vresid(1) is transferred to i_vresid(2) at output.
1012 : : !! The input f_atm(:,:,i_vresid(2)) contains minus the old HF forces.
1013 : : !! the value of i_vresid(2) is transferred to i_vresid(3) at output.
1014 : : !! The input f_atm(:,:,i_vresid(3)) contains minus the previous old HF forces.
1015 : : !! For the preconditioned change of atomic positions, the same logic as for the
1016 : : !! the potential residual is used, with i_vrespc replacing i_vresid.
1017 : : !! The input/output search vector is stored in f_atm(:,:,6)
1018 : : !! i_rhor(2:3)=index of the density (past and previous past) in the array f_fftgr
1019 : : !! i_vresid(3)=index of the residual potentials (present, past and previous
1020 : : !! past) in the array f_fftgr; also similar index for minus Hellman-Feynman
1021 : : !! forces in the array f_atm .
1022 : : !! i_vrespc(3)=index of the preconditioned residual potentials
1023 : : !! (present, past and previous past) in the array f_fftgr ;
1024 : : !! also similar index for the preconditioned change of atomic position (dtn_pc).
1025 : : !!
1026 : : !! TODO
1027 : : !! This routine is much too difficult to read ! Should be rewritten ...
1028 : : !! Maybe make separate subroutines for line search and CG step ?!
1029 : : !!
1030 : : !! SOURCE
1031 : :
1032 : 186 : subroutine scfcge(cplex,dbl_nnsclo,dtn_pc,etotal,f_atm,&
1033 : 186 : & f_fftgr,initialized,iscf,isecur,istep,&
1034 : 186 : & i_rhor,i_vresid,i_vrespc,moved_atm_inside,mpicomm,mpi_summarize,&
1035 : 186 : & natom,nfft,nfftot,nspden,n_fftgr,n_index,opt_denpot,response,rhor,ucvol,vtrial,xred,errid,errmess)
1036 : :
1037 : : !Arguments ------------------------------------
1038 : : !scalars
1039 : : integer,intent(in) :: cplex,initialized,iscf,isecur,istep,moved_atm_inside,mpicomm
1040 : : integer,intent(in) :: n_fftgr,n_index,natom,nfft,nfftot,nspden,opt_denpot,response
1041 : : integer,intent(out) :: dbl_nnsclo, errid
1042 : : character(len = 500), intent(out) :: errmess
1043 : : logical, intent(in) :: mpi_summarize
1044 : : real(dp),intent(in) :: etotal,ucvol
1045 : : !arrays
1046 : : integer,intent(inout) :: i_rhor(n_index),i_vresid(n_index),i_vrespc(n_index)
1047 : : real(dp),intent(in) :: dtn_pc(3,natom),rhor(cplex*nfft,nspden)
1048 : : real(dp),intent(inout) :: f_atm(3,natom,n_fftgr)
1049 : : real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr)
1050 : : real(dp),intent(inout) :: vtrial(cplex*nfft,nspden),xred(3,natom)
1051 : :
1052 : : !Local variables-------------------------------
1053 : : !mlinmin gives the maximum number of steps in the line minimization
1054 : : ! after which the algorithm is restarted (with a decrease of the
1055 : : ! adaptative trial step length). This number should not be large,
1056 : : ! since if the potential landscape is harmonic, the number of
1057 : : ! search steps should be small. If it is large, we are not in the
1058 : : ! harmonic region, and the CG algorithm will not be really useful,
1059 : : ! so one can just restart the algorithm ...
1060 : : !scalars
1061 : : integer,parameter :: mlinmin=5
1062 : : integer,save :: end_linmin,iline_cge,ilinear,ilinmin,isecur_eff,nlinear
1063 : : integer,save :: number_of_restart,status
1064 : : integer :: choice,iatom,idir,ifft,iline_cge_input,ilinmin_input,isp
1065 : : integer :: testcg,tmp,errid_
1066 : : real(dp),save :: d2edv2_old2,d_lambda_old2,dedv_old2,etotal_old
1067 : : real(dp),save :: etotal_previous=MAGIC_UNDEF,lambda_adapt,lambda_new,lambda_old,resid_old
1068 : : real(dp) :: d2e11,d2e12,d2e22,d2edv2_new,d2edv2_old
1069 : : real(dp) :: d2edv2_predict,d_lambda,de1,de2,dedv_mix
1070 : : real(dp) :: dedv_new,dedv_old,dedv_predict,determ,etotal_input
1071 : : real(dp) :: etotal_predict,gamma,lambda_input,lambda_predict2
1072 : : real(dp) :: lambda_predict=1.0_dp,ratio,reduction
1073 : : real(dp) :: resid_input,temp
1074 : : character(len=500) :: message
1075 : : !arrays
1076 : : real(dp) :: resid_new(1)
1077 : 186 : real(dp), allocatable :: tmp_fft1(:,:)
1078 : :
1079 : : ! *************************************************************************
1080 : :
1081 : 186 : errid = AB7_NO_ERROR
1082 : 186 : dbl_nnsclo = 0
1083 : :
1084 : : !reduction gives the level of reduction of the error in
1085 : : !the line minimization to be reached for the minimization to be
1086 : : !considered successfull
1087 : 186 : reduction=0.1_dp
1088 : :
1089 : : !nlinear increases with the number of times the 2D minimization succeded
1090 : : !to reach the true minimum directly. It is a measure of the
1091 : : !degree of parabolicity of the problem, and is used to
1092 : : !skip some steps by performing extrapolation.
1093 [ + + ]: 186 : if(istep==1)then
1094 : :
1095 : : ! Skipping some steps is sometimes unsecure, so it is possible
1096 : : ! to make nlinear start at a negative value - if isecur is positive
1097 : 30 : isecur_eff=isecur
1098 : 30 : nlinear=min(-isecur_eff,0)
1099 : 30 : ilinear=0
1100 : :
1101 : : ! Response function calculation are intrinsically harmonic, so one
1102 : : ! can shift isecur (by -2), and start with a positive nlinear
1103 [ - + ]: 30 : if(response==1)then
1104 : 0 : isecur_eff=isecur-2
1105 : 0 : nlinear=-isecur_eff
1106 : 0 : ilinear=nlinear
1107 : : end if
1108 : :
1109 : 30 : iline_cge=0
1110 : 30 : ilinmin=0
1111 : : end if
1112 : :
1113 : : !Compute actual residual resid_new (residual of f_fftgr(:,:,i_vrespc(1))
1114 : 186 : call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr)
1115 : :
1116 : : !Save input residual and ilinmin for final printing
1117 : 186 : resid_input=resid_new(1)
1118 : 186 : etotal_input=etotal
1119 : 186 : ilinmin_input=ilinmin
1120 : 186 : iline_cge_input=iline_cge
1121 : : !Transfer dtn_pc in f_atm
1122 [ + + ]: 186 : if(moved_atm_inside==1)then
1123 [ + + + + ]: 180 : f_atm(:,:,i_vrespc(1))=dtn_pc(:,:)
1124 : : end if
1125 : :
1126 : : !=======================================================================
1127 : : !Now the routine is decomposed in three mutually exclusive parts :
1128 : : !if(istep==1)then initialize the algorithm
1129 : : !else if(ilinmin>0)then perform the line minimisation
1130 : : !else if(ilinmin==0)then determine the new search direction (CG step)
1131 : : !=======================================================================
1132 : :
1133 : :
1134 : : !--------------------------------------
1135 : : !Here initialize the algorithm
1136 [ + + ]: 186 : if(istep==1)then
1137 : :
1138 : : ! At the beginning of each gstate run, lambda_adapt is forced to have the
1139 : : ! same value, that is 1.0_dp. In the other cases when istep=1 (at different
1140 : : ! broyden steps, for example), the previously obtained
1141 : : ! adaptive value is kept.
1142 [ + + ]: 30 : if(initialized==0)lambda_adapt=1.0_dp
1143 : 30 : lambda_old=0.0_dp
1144 : 30 : lambda_input=0.0_dp
1145 : 30 : number_of_restart=0
1146 : 30 : lambda_new=lambda_adapt
1147 : :
1148 [ + + + + ]: 310858 : f_fftgr(:,:,1)=vtrial(:,:)
1149 [ + + + + ]: 310858 : f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1150 : :
1151 : : ! This copy must be written in F77, because of stack problems on the DECs
1152 [ + + ]: 69 : do isp=1,nspden
1153 [ + + ]: 310858 : do ifft=1,cplex*nfft
1154 : 310828 : f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1))
1155 : : end do
1156 : : end do
1157 [ + + + + ]: 310858 : vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6)
1158 [ + + ]: 30 : if(moved_atm_inside==1)then
1159 [ + + + + ]: 27 : f_atm(:,:,1)=xred(:,:)
1160 [ + + + + ]: 27 : f_atm(:,:,i_rhor(2))=xred(:,:)
1161 : : ! There shouldn t be problems with the stack size for this small array.
1162 [ + + + + ]: 27 : f_atm(:,:,6)=f_atm(:,:,i_vrespc(1))
1163 [ + + + + ]: 27 : xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6)
1164 : : end if
1165 : 30 : tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
1166 : 30 : tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1167 : 30 : ilinmin=1
1168 : 30 : resid_old=resid_new(1)
1169 : 30 : etotal_old=etotal
1170 : :
1171 : 30 : status=0
1172 : :
1173 : : ! --------------------------------------
1174 : :
1175 : : ! Here performs the line minimisation
1176 [ + + ]: 156 : else if(ilinmin>0)then
1177 : :
1178 : 91 : lambda_input=lambda_new
1179 : :
1180 : : ! The choice with the Brent algorithm has been abandoned in version 1.6.m
1181 : :
1182 : : ! Compute the approximate energy derivatives dedv_new and dedv_old,
1183 : : ! from vresid and vresid_old
1184 : 91 : choice=2
1185 : : call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,&
1186 : : & f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,&
1187 : 91 : & natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred)
1188 : 91 : d_lambda=lambda_new-lambda_old
1189 : 91 : dedv_old=dedv_old/d_lambda
1190 : 91 : dedv_new=dedv_new/d_lambda
1191 : :
1192 : : ! DEBUG
1193 : : ! write(std_out,'(a,4es12.4,i3)' )' scfcge:lold,lnew,dold,dnew,status', &
1194 : : ! & lambda_old,lambda_new,dedv_old,dedv_new,status
1195 : : ! ENDDEBUG
1196 : :
1197 [ + - ]: 91 : if(status==0 .or. status==3)then
1198 : : !
1199 : : ! Then, compute a predicted point along the line
1200 : : ! The value of choice determines the minimization algorithm
1201 : : ! choice=1 uses the two values of the derivative of the energy
1202 : : ! choice=2 uses the two values of the energy, and and estimate of the
1203 : : ! second derivative at the mid-point.
1204 : :
1205 : 91 : choice=1
1206 [ + + ]: 91 : if(iscf==6)choice=2
1207 : : call findminscf(choice,dedv_new,dedv_old,dedv_predict,&
1208 : : & d2edv2_new,d2edv2_old,d2edv2_predict,&
1209 : : & etotal,etotal_old,etotal_predict,&
1210 : 91 : & lambda_new,lambda_old,lambda_predict,errid_,message)
1211 [ + + ]: 91 : if (errid_ /= AB7_NO_ERROR) then
1212 : 4 : call wrtout(std_out,message,'COLL')
1213 : : end if
1214 : :
1215 : : ! Suppress the next line for debugging (there is another such line)
1216 : 91 : status=0
1217 : :
1218 : : ! DEBUG
1219 : : ! Keep this debugging feature : it gives access to the investigation of lines
1220 : : ! in a different approach
1221 : : ! if(response==1 .and. istep>8)then
1222 : : ! lambda_predict=1.2d-2
1223 : : ! if(istep>=15)lambda_predict=lambda_predict-0.002
1224 : : ! if(istep>=14)stop
1225 : : ! status=3
1226 : : ! end if
1227 : : ! ENDDEBUG
1228 : :
1229 : : else
1230 [ # # ]: 0 : if(status/=-1)then
1231 : 0 : status=-1
1232 : 0 : lambda_predict=-2.5_dp
1233 : : else
1234 : 0 : lambda_predict=lambda_predict+0.1_dp
1235 : : end if
1236 : : end if
1237 : :
1238 : : ! If the predicted point is very close to the most recent
1239 : : ! computed point, while this is the first trial on this line,
1240 : : ! then we are in the linear regime :
1241 : : ! nlinear is increased by one unit. For the time being, do this even when
1242 : : ! moved_atm_inside==1 (the code still works when it is done, but it
1243 : : ! seems to be a bit unstable). The maximal value of nlinear is 1, except
1244 : : ! when isecur_eff is a negative number, less than -1.
1245 : : if( abs(lambda_predict-lambda_new)/&
1246 [ + + + + ]: 91 : & (abs(lambda_predict)+abs(lambda_new)) < 0.01 .and. ilinmin==1 ) then
1247 : : ! if(moved_atm_inside==0 .and. nlinear<max(1,-isecur_eff) )nlinear=nlinear+1
1248 [ + + ]: 17 : if(nlinear<max(1,-isecur_eff))nlinear=nlinear+1
1249 : 17 : ilinear=nlinear
1250 : : end if
1251 : :
1252 : : ! If the predicted point is close to the most recent computed point,
1253 : : ! or the previous one, set on the flag of end of line minization
1254 : 91 : end_linmin=0
1255 [ + + ]: 91 : if(abs(lambda_new-lambda_predict)*2.0_dp&
1256 : 54 : & /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1
1257 [ + + ]: 91 : if(abs(lambda_old-lambda_predict)*2.0_dp&
1258 : 3 : & /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1
1259 : :
1260 [ - + ]: 91 : if(status/=0)end_linmin=0
1261 : :
1262 : : ! Save the closest old lambda, if needed,
1263 : : ! also examine the reduction of the interval, and eventual stop
1264 : : ! the present line minimisation, because of convergence (end_linmin=1)
1265 : : ! Also treat the case in which the predicted value of lambda is negative,
1266 : : ! or definitely too small in which case the algorithm has to be restarted
1267 : : ! (not a very good solution, though ...)
1268 : : ! Finally also treat the case where insufficiently converged
1269 : : ! density at lambda=0.0_dp happens, which screws up the line minimisation.
1270 : :
1271 : : ! Here restart the algorithm with the best vtrial.
1272 : : ! Also make reduction in lambda_adapt
1273 : : ! DEBUG
1274 : : ! write(std_out,*)' scfcge : status=',status
1275 : : ! ENDDEBUG
1276 [ + + + - : 91 : if( end_linmin==0 .and. status==0 .and. &
- + - - -
+ - - -
+ ]
1277 : : & ( (lambda_predict<0.005_dp*lambda_adapt .and. iscf==5) .or. &
1278 : : & (abs(lambda_predict)<0.005_dp*lambda_adapt .and. iscf==6).or. &
1279 : : & ilinmin==mlinmin ) )then
1280 [ # # ]: 0 : if(number_of_restart>12)then
1281 : 0 : errid = AB7_ERROR_MIXING_CONVERGENCE
1282 : : write(errmess,'(a,a,i0,a,a,a,a,a)')&
1283 : 0 : & 'Potential-based CG line minimization not',' converged after ',number_of_restart,' restarts. ',ch10,&
1284 : 0 : & 'Action : read the eventual warnings about lack of convergence.',ch10,&
1285 : 0 : & 'Some might be relevant. Otherwise, raise nband. Returning'
1286 : 0 : ABI_WARNING(errmess)
1287 : 0 : return
1288 : : end if
1289 : : ! Make reduction in lambda_adapt (kind of steepest descent...)
1290 : : write(message,'(a,a,a)')&
1291 : 0 : & 'Potential-based CG line minimization has trouble to converge.',ch10,&
1292 : 0 : & 'The algorithm is restarted with more secure parameters.'
1293 : 0 : ABI_WARNING(message)
1294 : 0 : number_of_restart=number_of_restart+1
1295 : : ! At the second restart, double the number of non-self consistent loops.
1296 [ # # ]: 0 : if(number_of_restart>=2)dbl_nnsclo=1
1297 : 0 : lambda_adapt=lambda_adapt*0.7_dp
1298 : 0 : lambda_new=lambda_adapt
1299 : : ! If the last energy is better than the old one, transfer the data.
1300 : : ! Otherwise, no transfer must occur (very simple to code...)
1301 [ # # # # ]: 0 : if(etotal<etotal_old .or. abs(lambda_old)<1.0d-8)then
1302 [ # # # # ]: 0 : f_fftgr(:,:,1)=vtrial(:,:)
1303 [ # # # # ]: 0 : f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1304 [ # # ]: 0 : do isp=1,nspden
1305 [ # # ]: 0 : do ifft=1,cplex*nfft
1306 : 0 : f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1))
1307 : : end do
1308 : : end do
1309 [ # # ]: 0 : if(moved_atm_inside==1)then
1310 [ # # # # ]: 0 : f_atm(:,:,1)=xred(:,:)
1311 [ # # # # ]: 0 : f_atm(:,:,i_rhor(2))=xred(:,:)
1312 [ # # # # ]: 0 : f_atm(:,:,6)=f_atm(:,:,i_vrespc(1))
1313 : : end if
1314 : 0 : tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
1315 : 0 : tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1316 : 0 : resid_old=resid_new(1)
1317 : 0 : etotal_old=etotal
1318 : : end if
1319 : 0 : lambda_old=0.0_dp
1320 : 0 : ilinmin=1
1321 : : ! Putting the flag to -1 avoids the usual actions taken with end_linmin=1
1322 : 0 : end_linmin=-1
1323 : : ! Also put ilinear and nlinear to 0
1324 : 0 : ilinear=0
1325 : 0 : nlinear=0
1326 : :
1327 : : ! Here lambda_new is the closest to lambda_predict,
1328 : : ! or lambda_old is still 0.0_dp, while the energy shows that the minimum
1329 : : ! is away from 0.0_dp (insufficiently converged density at lambda=0.0_dp).
1330 : : else if( abs(lambda_new-lambda_predict)<abs(lambda_old-lambda_predict) &
1331 [ + + + + : 91 : & .or. &
- + - - ]
1332 : : & ( abs(lambda_old)<1.0d-6 .and. &
1333 : : & ilinmin>1 .and. &
1334 : : & etotal>etotal_previous ) &
1335 : : & )then
1336 [ + + + + ]: 950947 : f_fftgr(:,:,1)=vtrial(:,:)
1337 : 84 : tmp=i_rhor(3) ; i_rhor(3)=i_rhor(2) ; i_rhor(2)=tmp
1338 [ + + + + ]: 950947 : f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1339 : 84 : tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2)
1340 : 84 : i_vrespc(2)=i_vrespc(1); i_vrespc(1)=tmp;
1341 : 84 : tmp=i_vresid(3); i_vresid(3)=i_vresid(2)
1342 : 84 : i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1343 [ + + ]: 84 : if(moved_atm_inside==1)then
1344 [ + + + + ]: 90 : f_atm(:,:,1)=xred(:,:)
1345 [ + + + + ]: 90 : f_atm(:,:,i_rhor(2))=xred(:,:)
1346 : : end if
1347 : 84 : d_lambda_old2=lambda_old-lambda_new
1348 : 84 : lambda_old=lambda_new
1349 : 84 : etotal_old=etotal
1350 : 84 : resid_old=resid_new(1)
1351 : 84 : d2edv2_old2=d2edv2_new
1352 : 84 : dedv_old=dedv_new
1353 : 84 : dedv_old2=dedv_new
1354 : : ! if(abs(lambda_new-lambda_predict)*2.0_dp&
1355 : : ! & /abs(lambda_new+lambda_predict) <reduction) end_linmin=1
1356 : :
1357 : : ! Here lambda_old is the closest to lambda_predict (except for avoiding
1358 : : ! lambda_old==0.0_dp)
1359 : : else
1360 : 7 : tmp=i_vresid(3) ; i_vresid(3)=i_vresid(1) ; i_vresid(1)=tmp
1361 [ + + + + ]: 76017 : f_fftgr(:,:,i_rhor(3))=rhor(:,:)
1362 [ + + + + : 15 : if(moved_atm_inside==1) f_atm(:,:,i_rhor(3))=xred(:,:)
+ + ]
1363 : 7 : tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(1) ; i_vrespc(1)=tmp
1364 : 7 : d_lambda_old2=lambda_new-lambda_old
1365 : 7 : etotal_previous=etotal
1366 : 7 : d2edv2_old2=d2edv2_old
1367 : 7 : dedv_old2=dedv_old
1368 : : ! if(abs(lambda_old-lambda_predict)*2.0_dp&
1369 : : ! & /abs(lambda_old+lambda_predict) <reduction) end_linmin=1
1370 : : end if
1371 : :
1372 : : ! If the interval has not yet been sufficiently reduced,
1373 : : ! continue the search
1374 [ + + ]: 91 : if(end_linmin==0)then
1375 : 34 : lambda_new=lambda_predict
1376 : :
1377 : : ! DEBUG
1378 : : ! write(std_out,'(a,2es16.6)' )&
1379 : : ! & ' scfcge : continue search, lambda_old,lambda_new=',lambda_old,lambda_new
1380 : : ! write(std_out,'(a,2es16.6)' )&
1381 : : ! & ' scfcge : f_fftgr(3:4,1,1)=',f_fftgr(3:4,1,1)
1382 : : ! write(std_out,'(a,2es16.6)' )&
1383 : : ! & ' scfcge : f_fftgr(3:4,1,6)=',f_fftgr(3:4,1,6)
1384 : : ! ENDDEBUG
1385 : :
1386 [ + + + + ]: 327307 : vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6)
1387 [ + + ]: 34 : if(moved_atm_inside==1)then
1388 [ + + + + ]: 45 : xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6)
1389 : : end if
1390 : :
1391 : 34 : ilinmin=ilinmin+1
1392 : : !
1393 : : ! Here generates a starting point for next line search
1394 : : else
1395 : 57 : iline_cge=iline_cge+1
1396 [ + - ]: 57 : if(end_linmin==1)ilinmin=0
1397 : 57 : lambda_old=0.0_dp
1398 : :
1399 : : ! In order to generate the new step, take into account previous
1400 : : ! optimal lambdas (including those of previous ion moves),
1401 : : ! and the selected new one, if it is positive.
1402 : : ! However, wait iline_cge>1 to select new ones.
1403 : : ! lambda_adapt has been initialized at 1.0_dp
1404 [ + + + - ]: 57 : if(iline_cge>1 .and. lambda_new>0.0_dp )then
1405 : : ! Actually compute a geometric mean
1406 : : lambda_adapt= ( lambda_adapt**(dble(iline_cge-1)) * abs(lambda_new)) &
1407 : 29 : & **(1.0_dp/dble(iline_cge))
1408 : : ! In order to recover the previous algorithm, it is enough
1409 : : ! to decomment the next line
1410 : : ! lambda_adapt=1.0_dp
1411 : : end if
1412 : 57 : lambda_new=lambda_adapt
1413 : :
1414 [ + + + + ]: 699657 : vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2))
1415 [ + + ]: 57 : if(moved_atm_inside==1)then
1416 [ + + + + ]: 54 : xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2))
1417 : : end if
1418 : :
1419 : : ! End choice between continue line minim and determine new direction
1420 : : end if
1421 : :
1422 : : !
1423 : : ! -------------------------------
1424 : :
1425 : : ! Here perform the CG step
1426 : :
1427 [ + - ]: 65 : else if(ilinmin==0)then
1428 : :
1429 : : ! Compute the approximate energy derivatives dedv_mix,dedv_new,dedv_old
1430 : 65 : choice=3
1431 : : call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,&
1432 : : & f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,&
1433 : 65 : & natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred)
1434 : :
1435 : 65 : dedv_mix=dedv_mix/lambda_new
1436 : 65 : dedv_new=dedv_new/lambda_new
1437 : 65 : dedv_old=dedv_old/lambda_new
1438 : :
1439 : : ! DEBUG
1440 : : ! write(message, '(a,3es12.4)' )' scfcge: lambda_adapt',&
1441 : : ! & lambda_adapt
1442 : : ! call wrtout(std_out,message,'COLL')
1443 : :
1444 : : ! write(message, '(a,3es12.4)' )' scfcge: dedv_old,dedv_new,dedv_mix',&
1445 : : ! & dedv_old,dedv_new,dedv_mix
1446 : : ! call wrtout(std_out,message,'COLL')
1447 : : ! ENDDEBUG
1448 : :
1449 : : ! Then, compute a predicted point, either along the line,
1450 : : ! or in a 2D plane
1451 : 65 : testcg=1
1452 : : if(testcg==0)then
1453 : : ! This part corresponds to steepest descent,
1454 : : ! in which the line minimisation can be done
1455 : : ! using different algorithms, varying with the value of choice
1456 : : choice=1
1457 : : if(iscf==6)choice=2
1458 : : call findminscf(choice,dedv_new,dedv_old,dedv_predict,&
1459 : : & d2edv2_new,d2edv2_old,d2edv2_predict,&
1460 : : & etotal,etotal_old,etotal_predict,&
1461 : : & lambda_new,lambda_old,lambda_predict,errid_,message)
1462 : : if (errid_ /= AB7_NO_ERROR) then
1463 : : call wrtout(std_out,message,'COLL')
1464 : : end if
1465 : : lambda_predict2=0.0_dp
1466 : : ! Suppress the next line for debugging (there is another such line)
1467 : : status=0
1468 : : else
1469 : : ! This part corresponds to conjugate gradient
1470 : : ! A 2D minimisation is performed
1471 : : ! oldest direction is labelled 2
1472 : : ! newest direction is labelled 1
1473 : 65 : de1=dedv_old ; de2=dedv_old2
1474 : 65 : d2e11=(dedv_new-dedv_old)/lambda_new
1475 : 65 : d2e22=d2edv2_old2
1476 : 65 : d2e12=(dedv_mix-dedv_old)/d_lambda_old2
1477 : : ! The system to be solved is
1478 : : ! 0 = de1 + lambda1 d2e11 + lambda2 d2d12
1479 : : ! 0 = de2 + lambda1 d2e12 + lambda2 d2d22
1480 : 65 : determ=d2e11*d2e22-d2e12*d2e12
1481 : 65 : lambda_predict=-(de1*d2e22-de2*d2e12)/determ
1482 : 65 : lambda_predict2=(de1*d2e12-de2*d2e11)/determ
1483 : 65 : d2edv2_new=d2e11 ; d2edv2_old=d2e11
1484 : : end if
1485 : :
1486 : : ! DEBUG
1487 : : ! write(message, '(a,5es11.3)' )' scfcge: de1,de2,d2e11,d2e22,d2e12',&
1488 : : ! & de1,de2,d2e11,d2e22,d2e12
1489 : : ! call wrtout(std_out,message,'COLL')
1490 : : ! write(std_out,'(a,2es12.4)' )' scfcge: la_predict,la_predict2',&
1491 : : ! & lambda_predict,lambda_predict2
1492 : : ! -----
1493 : : ! write(std_out,*)'residues ',
1494 : : ! !$ de1+lambda_predict*d2e11+lambda_predict2*d2e12,
1495 : : ! !$ de2+lambda_predict*d2e12+lambda_predict2*d2e22
1496 : : ! if(.true.)stop
1497 : : ! ENDDEBUG
1498 : : !
1499 : :
1500 : : ! Determine the region of the 2D search space
1501 : : ! in which the predicted point is located,
1502 : : ! or use linear indicator to decide interpolation
1503 : : ! and advance to next 2D search.
1504 : 65 : end_linmin=0
1505 : 65 : write(message, '(a,2i3)' )' nlinear, ilinear',nlinear,ilinear
1506 : 65 : call wrtout(std_out,message,'COLL')
1507 [ - + ]: 65 : if(lambda_predict<0.0_dp)then
1508 : : ! Something is going wrong. Just take a reasonable step
1509 : : ! along the steepest descent direction (Region III).
1510 : : ! Actually, Region I and region III are treated in the same way later.
1511 : : ! In effect, this corresponds to restart the algorithm
1512 : 0 : end_linmin=3
1513 : : ! Also put ilinear and nlinear to 0
1514 : 0 : ilinear=0
1515 : 0 : nlinear=0
1516 : : ! Decrease the adaptive step to predict next direction
1517 : 0 : lambda_adapt=lambda_adapt*0.7_dp
1518 [ + + ]: 65 : else if(ilinear>=1) then
1519 : : ! Region IV : will do an interpolation
1520 : 6 : end_linmin=4
1521 : 6 : ilinear=ilinear-1
1522 : : else if(abs(lambda_predict2)>reduction .or.&
1523 : : & lambda_predict<0.5_dp .or.&
1524 [ + - + + : 42 : & lambda_predict>2.5_dp .or.&
- + ]
1525 [ + + ]: 59 : & lambda_predict-abs(lambda_predict2)/reduction <0.0_dp ) then
1526 : : ! Region II : lambda_predict is not too good, and not too bad.
1527 : 18 : end_linmin=2
1528 [ + + ]: 41 : else if (abs(1.0_dp-lambda_predict)<reduction)then
1529 : : ! Region I, the out-of-line point is OK.
1530 : 13 : end_linmin=1
1531 : : else
1532 : : ! If everything fails, then region II.
1533 : 28 : end_linmin=2
1534 : : end if
1535 : :
1536 : : ! DEBUG
1537 : : ! write(message, '(a,2es12.4,i2)' )&
1538 : : ! & ' scfcge : la_predict, la_predict2, region',&
1539 : : ! & lambda_predict,lambda_predict2,end_linmin
1540 : : ! call wrtout(std_out,message,'COLL')
1541 : : ! ENDDEBUG
1542 : :
1543 : : ! Treat region I, in the same way as region III
1544 [ + + ]: 65 : if(end_linmin==1 .or. end_linmin==3)then
1545 : :
1546 : : ! In region I, the line search is
1547 : : ! along vtrial-vtrial_old.
1548 : : ! The closest point is the new point
1549 : : ! thus to be transfered in the "old" locations
1550 : :
1551 [ + + ]: 26 : do isp=1,nspden
1552 [ + + ]: 53568 : do ifft=1,cplex*nfft
1553 : 53555 : f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new
1554 : : end do
1555 : : end do
1556 [ + + + + ]: 53568 : f_fftgr(:,:,1)=vtrial(:,:)
1557 [ + + + + ]: 53568 : f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1558 [ - + ]: 13 : if(moved_atm_inside==1)then
1559 [ # # # # ]: 0 : f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new
1560 [ # # # # ]: 0 : f_atm(:,:,1)=xred(:,:)
1561 [ # # # # ]: 0 : f_atm(:,:,i_rhor(2))=xred(:,:)
1562 : : end if
1563 : 13 : tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
1564 : 13 : tmp=i_vresid(3) ; i_vresid(3)=i_vresid(2)
1565 : 13 : i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1566 : 13 : d_lambda_old2=-lambda_new
1567 : : lambda_old=lambda_new
1568 : 13 : etotal_old=etotal
1569 : 13 : resid_old=resid_new(1)
1570 : : d2edv2_old=d2edv2_new
1571 : 13 : dedv_old=dedv_new
1572 : :
1573 : : ! Region I or III : one is close of the 2D minimum,
1574 : : ! or lambda_predict was negative (indicate a problem of convergence)
1575 : : ! Compute next trial potential along the
1576 : : ! PC residual and not along this search direction.
1577 : 13 : ilinmin=0
1578 : : ! Question : isn t it here that one should prevent region I to called
1579 : : ! itself more than 1 time ???
1580 : : ! Here the small difference between region I and region III
1581 [ - + ]: 13 : if(end_linmin==3)ilinmin=1
1582 : 13 : lambda_old=0.0_dp
1583 : 13 : lambda_new=lambda_adapt
1584 : :
1585 [ + + + + ]: 53568 : vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2))
1586 [ - + ]: 13 : if(moved_atm_inside==1)then
1587 [ # # # # ]: 0 : xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2))
1588 : : end if
1589 : : ! The new vtrial has been generated
1590 : :
1591 : : else
1592 : :
1593 : : ! Here region II or IV
1594 : 52 : ilinmin=1
1595 [ + - ]: 52 : if (lambda_predict==0._dp) then
1596 : : gamma=zero
1597 : : else
1598 : 52 : gamma=lambda_predict2/lambda_predict
1599 : : end if
1600 : : ! Compute new search direction and trial potential
1601 : 52 : write(message,*)' compute new search direction '
1602 : 52 : call wrtout(std_out,message,'COLL')
1603 [ + + ]: 124 : do isp=1,nspden
1604 [ + + ]: 676792 : do ifft=1,cplex*nfft
1605 : : f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new+ &
1606 : 676740 : & gamma*f_fftgr(ifft,isp,6)
1607 : : end do
1608 : : end do
1609 [ + + + + ]: 676792 : vtrial(:,:)=f_fftgr(:,:,1)+ lambda_predict*f_fftgr(:,:,6)
1610 [ + + ]: 52 : if(moved_atm_inside==1)then
1611 [ + + + + ]: 54 : f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new+ gamma*f_atm(:,:,6)
1612 [ + + + + ]: 54 : xred(:,:)=f_atm(:,:,1)+ lambda_predict*f_atm(:,:,6)
1613 : : end if
1614 : :
1615 : : ! If end_linmin==2, then this vtrial is the good one
1616 : :
1617 [ + + ]: 52 : if(end_linmin==2)then
1618 : :
1619 : 46 : lambda_old=0.0_dp
1620 : 46 : lambda_new=lambda_predict
1621 : :
1622 [ + - ]: 6 : else if(end_linmin==4)then
1623 : :
1624 : : ! predict the result of the computation at the trial potential
1625 : : ! defined in the end_linmin==2 case
1626 : 6 : gamma=lambda_predict2/d_lambda_old2
1627 : 6 : ratio=lambda_predict/lambda_new
1628 : :
1629 : : ! Take care of vtrial
1630 [ + + + + ]: 84679 : f_fftgr(:,:,1)=vtrial(:,:)
1631 : :
1632 [ - + + - : 24 : ABI_MALLOC(tmp_fft1,(cplex*nfft,nspden))
- + + - +
- + - - +
- + ]
1633 : : ! Take care of vresid
1634 [ + + + + ]: 84679 : tmp_fft1(:,:)=f_fftgr(:,:,i_vresid(2))
1635 : : f_fftgr(:,:,i_vresid(2))=tmp_fft1(:,:)&
1636 : : & +ratio*(f_fftgr(:,:,i_vresid(1))-tmp_fft1(:,:))&
1637 [ + + + + ]: 84679 : & +gamma*(f_fftgr(:,:,i_vresid(3))-tmp_fft1(:,:))
1638 [ + + + + ]: 84679 : f_fftgr(:,:,i_vresid(3))=tmp_fft1(:,:)
1639 : :
1640 : : ! Take care of rhor
1641 [ + + + + ]: 84679 : tmp_fft1(:,:)=f_fftgr(:,:,i_rhor(2))
1642 : : f_fftgr(:,:,i_rhor(2))=tmp_fft1(:,:)&
1643 : : & +ratio*(rhor(:,:)-tmp_fft1(:,:))&
1644 [ + + + + ]: 84679 : & +gamma*(f_fftgr(:,:,i_rhor(3))-tmp_fft1(:,:))
1645 [ + + + + ]: 84679 : f_fftgr(:,:,i_rhor(3))=tmp_fft1(:,:)
1646 : :
1647 : : ! Take care of vrespc
1648 [ + + + + ]: 84679 : tmp_fft1(:,:)=f_fftgr(:,:,i_vrespc(2))
1649 : : f_fftgr(:,:,i_vrespc(2))=tmp_fft1(:,:)&
1650 : : & +ratio*(f_fftgr(:,:,i_vrespc(1))-tmp_fft1(:,:))&
1651 [ + + + + ]: 84679 : & +gamma*(f_fftgr(:,:,i_vrespc(3))-tmp_fft1(:,:))
1652 [ + + + + ]: 84679 : f_fftgr(:,:,i_vrespc(3))=tmp_fft1(:,:)
1653 [ - + ]: 6 : ABI_FREE(tmp_fft1)
1654 : :
1655 [ + + ]: 6 : if(moved_atm_inside==1)then
1656 [ + + ]: 8 : do idir=1,3
1657 [ + + ]: 20 : do iatom=1,natom
1658 : :
1659 : : ! Take care of xred
1660 : 12 : f_atm(idir,iatom,1)=xred(idir,iatom)
1661 : :
1662 : : ! Take care of -HF forces
1663 : 12 : temp=f_atm(idir,iatom,i_vresid(2))
1664 : : f_atm(idir,iatom,i_vresid(2))=f_atm(idir,iatom,i_vresid(2))&
1665 : : & +ratio*(f_atm(idir,iatom,i_vresid(1))-f_atm(idir,iatom,i_vresid(2)))&
1666 : 12 : & +gamma*(f_atm(idir,iatom,i_vresid(3))-f_atm(idir,iatom,i_vresid(2)))
1667 : 12 : f_atm(idir,iatom,i_vresid(3))=temp
1668 : :
1669 : : ! Take care of old xreds
1670 : 12 : temp=f_atm(idir,iatom,i_rhor(2))
1671 : : f_atm(idir,iatom,i_rhor(2))=f_atm(idir,iatom,i_rhor(2))&
1672 : : & +ratio*( xred(idir,iatom) -f_atm(idir,iatom,i_rhor(2)))&
1673 : 12 : & +gamma*(f_atm(idir,iatom,i_rhor(3))-f_atm(idir,iatom,i_rhor(2)))
1674 : 12 : f_atm(idir,iatom,i_rhor(3))=temp
1675 : :
1676 : : ! Take care of preconditioned changes of atomic positions
1677 : 12 : temp=f_atm(idir,iatom,i_vrespc(2))
1678 : : f_atm(idir,iatom,i_vrespc(2))=f_atm(idir,iatom,i_vrespc(2))&
1679 : : & +ratio*(f_atm(idir,iatom,i_vrespc(1))-f_atm(idir,iatom,i_vrespc(2)))&
1680 : 12 : & +gamma*(f_atm(idir,iatom,i_vrespc(3))-f_atm(idir,iatom,i_vrespc(2)))
1681 : 18 : f_atm(idir,iatom,i_vrespc(3))=temp
1682 : :
1683 : : end do
1684 : : end do
1685 : : end if
1686 : :
1687 : : ! Since we are at the 2D minimum, the derivative is supposed
1688 : : ! to vanish. Note that dedv_old should not change, by contrast.
1689 : 6 : dedv_old2=0.0_dp
1690 : 6 : d_lambda_old2=-lambda_predict
1691 : 6 : d2edv2_old2=-dedv_old/lambda_predict
1692 : 6 : lambda_old=lambda_predict
1693 : 6 : ilinmin=0
1694 : :
1695 : : ! So, jump to the next line
1696 : 6 : iline_cge=iline_cge+1
1697 : 6 : write(message,*)' energy CG update : after 2D interpolation,'
1698 : 6 : call wrtout(std_out,message,'COLL')
1699 : 6 : write(message,*)' computation in the next plane '
1700 : 6 : call wrtout(std_out,message,'COLL')
1701 : 6 : write(message,*)
1702 : 6 : call wrtout(std_out,message,'COLL')
1703 : 6 : lambda_old=0.0_dp
1704 : 6 : lambda_new=lambda_adapt
1705 : :
1706 [ + + + + ]: 84679 : vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2))
1707 [ + + ]: 6 : if(moved_atm_inside==1)then
1708 [ + + + + ]: 18 : xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2))
1709 : : end if
1710 : :
1711 : : ! The new trial potential is now generated
1712 : :
1713 : : ! End the specific treatment of region IV
1714 : : end if
1715 : : !
1716 : : ! End the choice between treatment of region I, II, or IV
1717 : : end if
1718 : :
1719 : : ! End of choice between initialisation or more developed parts of the CG algorithm
1720 : : else
1721 : 0 : errid = AB7_ERROR_MIXING_ARG
1722 : 0 : errmess = 'scfcge : BUG You should not be here ! '
1723 : 0 : return
1724 : : end if
1725 : :
1726 : : !--------------------------------------
1727 : :
1728 : : !Write information : it will be easy to read by typing grep scfcge logfile
1729 : :
1730 [ + + ]: 186 : if(istep==1)then
1731 : 30 : write(message,'(a,a,a)') ' scfcge:',ch10,' scfcge:istep-iline_cge-ilinmin lambda etot resid '
1732 : 30 : call wrtout(std_out,message,'COLL')
1733 : : end if
1734 : :
1735 [ + + + + ]: 186 : if(ilinmin_input/=0 .or. istep==1)then
1736 : : ! Usual line minimisation step
1737 : :
1738 [ + - ]: 121 : if(iline_cge_input<10)then
1739 : : write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' )&
1740 : 121 : & ' scfcge: actual ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input
1741 : : else
1742 : : write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' )&
1743 : 0 : & ' scfcge: actual ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input
1744 : : end if
1745 : 121 : call wrtout(std_out,message,'COLL')
1746 : :
1747 [ + + + + ]: 121 : if( (end_linmin==1.or.end_linmin==-1) .and. istep/=1 )then
1748 : :
1749 [ + - ]: 57 : if(end_linmin==1)then
1750 : : write(message, '(a,es13.4,a,i2,a,a)' )&
1751 : 57 : & ' scfcge: predict ',lambda_predict,&
1752 : 57 : & ' suff. close => next line, ilinear=',ilinear,ch10,&
1753 : 114 : & ' scfcge:'
1754 [ # # ]: 0 : else if(end_linmin==-1)then
1755 : : write(message, '(a,es13.4,a,a,a)' )&
1756 : 0 : & ' scfcge: predict ',lambda_predict,&
1757 : 0 : & ' restart the algorithm ',ch10,&
1758 : 0 : & ' scfcge:'
1759 : : end if
1760 : 57 : call wrtout(std_out,message,'COLL')
1761 : :
1762 [ + - ]: 57 : if(iline_cge_input<9)then
1763 : : write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' ) &
1764 : 57 : & ' scfcge: start ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old
1765 : : else
1766 : : write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' ) &
1767 : 0 : & ' scfcge: start ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old
1768 : : end if
1769 : 57 : call wrtout(std_out,message,'COLL')
1770 : :
1771 [ + + ]: 64 : else if(istep/=1) then
1772 : : write(message, '(a,es13.4,a)' )&
1773 : 34 : & ' scfcge: predict ',lambda_predict,&
1774 : 68 : & ' not close enough => continue minim.'
1775 : 34 : call wrtout(std_out,message,'COLL')
1776 : : end if
1777 : :
1778 : : else
1779 : : ! CG prediction
1780 [ + - ]: 65 : if(iline_cge_input<10)then
1781 : : write(message, '(a,i4,a,i1,a,es11.4,es20.12,es12.4,a,i1)' )&
1782 : 65 : & ' scfcge: actual ',istep,'-',iline_cge_input,'-off',&
1783 : 130 : & lambda_adapt,etotal_input,resid_input,', end=',end_linmin
1784 : : else
1785 : : write(message, '(a,i3,a,i2,a,es11.4,es20.12,es12.4,a,i1)' )&
1786 : 0 : & ' scfcge: actual ',istep,'-',iline_cge_input,'-off',&
1787 : 0 : & lambda_adapt,etotal_input,resid_input,', end=',end_linmin
1788 : : end if
1789 : 65 : call wrtout(std_out,message,'COLL')
1790 : :
1791 [ + + ]: 65 : if(end_linmin==4)then
1792 : 6 : write(message, '(a)' ) ' scfcge:'
1793 : 6 : call wrtout(std_out,message,'COLL')
1794 : : end if
1795 : :
1796 : : end if
1797 : :
1798 : : end subroutine scfcge
1799 : : !!***
1800 : :
1801 : : !!****f* ABINIT/scfeig
1802 : : !! NAME
1803 : : !! scfeig
1804 : : !!
1805 : : !! FUNCTION
1806 : : !! Compute the largest eigenvalue and eigenvector of the SCF cycle.
1807 : : !! A brute force algorithm is presently used.
1808 : : !!
1809 : : !! INPUTS
1810 : : !! istep= number of the step in the SCF cycle
1811 : : !! nfft=(effective) number of FFT grid points (for this processor)
1812 : : !! nspden=number of spin-density components
1813 : : !!
1814 : : !! OUTPUT
1815 : : !! (see side effects)
1816 : : !!
1817 : : !! SIDE EFFECTS
1818 : : !! vtrial0(nfft,nspden)= contains vtrial at istep == 1
1819 : : !! vtrial(nfft,nspden)= at input, it is the trial potential that gave vresid .
1820 : : !! at output, it is an updated trial potential
1821 : : !! vrespc(nfft,nspden)=the input preconditioned residual potential
1822 : : !! work(nfft,nspden,2)=work space
1823 : : !!
1824 : : !! SOURCE
1825 : :
1826 : 47 : subroutine scfeig(istep,nfft,nspden,vrespc,vtrial,vtrial0,work,errid,errmess)
1827 : :
1828 : : !Arguments ------------------------------------
1829 : : !scalars
1830 : : integer,intent(in) :: istep,nfft,nspden
1831 : : integer,intent(out) :: errid
1832 : : character(len = 500), intent(out) :: errmess
1833 : : !arrays
1834 : : real(dp),intent(inout) :: vtrial0(nfft,nspden),work(nfft,nspden,2)
1835 : : real(dp),intent(inout) :: vrespc(nfft,nspden)
1836 : : real(dp), intent(inout) :: vtrial(nfft,nspden)
1837 : :
1838 : : !Local variables-------------------------------
1839 : : !scalars
1840 : : integer :: ifft,isp
1841 : : real(dp) :: eigen_scf,factor,fix_resid,resid_new,resid_old
1842 : : character(len=500) :: message
1843 : :
1844 : : ! *************************************************************************
1845 : :
1846 : 47 : errid = AB7_NO_ERROR
1847 : :
1848 [ - + ]: 47 : if(nspden==4)then
1849 : 0 : errid = AB7_ERROR_MIXING_ARG
1850 : 0 : write(errmess, *) ' scfeig: does not work yet for nspden=4'
1851 : 0 : return
1852 : : end if
1853 : :
1854 : : !Set a fixed residual square for normalization of eigenvectors
1855 : 47 : fix_resid=1.0d-4
1856 : :
1857 : : !A few initialisations for the first istep
1858 [ + + ]: 47 : if(istep==1)then
1859 : :
1860 : : write(message, '(a,es12.4,a,a,a,a,a,a,a)' )&
1861 : 3 : & ' scfeig: fixed PC_residual square =',fix_resid,ch10,&
1862 : 3 : & ' Note that fixed resid should always be much larger',ch10,&
1863 : 3 : & ' than initial PC resid square, still sufficiently',ch10,&
1864 : 6 : & ' small to reduce anharmonic effects ',ch10
1865 : 3 : call wrtout(std_out,message,'COLL')
1866 : :
1867 : : ! Compute the preconditioned residual
1868 : 3 : resid_old=0.0_dp
1869 [ + + ]: 7 : do isp=1,nspden
1870 [ + + ]: 36631 : do ifft=1,nfft
1871 : 36628 : resid_old=resid_old+vrespc(ifft,isp)**2
1872 : : end do
1873 : : end do
1874 : 3 : write(message, '(a,es12.4)' )' scfeig: initial PC_residual square =',resid_old
1875 : 3 : call wrtout(std_out,message,'COLL')
1876 [ - + ]: 3 : if(resid_old>1.0d-8)then
1877 : 0 : errid = AB7_ERROR_MIXING_ARG
1878 : 0 : write(errmess,'(a,a,a,a,a,a,a,a,a,a)') ch10,&
1879 : 0 : & ' scfeig : ERROR -',ch10,&
1880 : 0 : & ' This value is not good enough to allow',ch10,&
1881 : 0 : & ' the computation of the eigenvectors of the SCF cycle.',ch10,&
1882 : 0 : & ' It should be better than 1.0d-8 .',ch10,&
1883 : 0 : & ' Action : improve the accuracy of your starting wavefunctions.'
1884 : 0 : return
1885 : : end if
1886 : :
1887 : : ! Also transfer vtrial in vtrial_old
1888 [ + + + + ]: 36631 : vtrial0(:,:)=vtrial(:,:)
1889 : :
1890 : : ! In order to start the search for eigenvectors,
1891 : : ! use the tiny residual vector, renormalized
1892 : 3 : factor=sqrt(fix_resid/resid_old)
1893 [ + + + + ]: 36631 : work(:,:,1)=vrespc(:,:)*factor
1894 [ + + + + ]: 36631 : vtrial(:,:)=vtrial0(:,:)+work(:,:,1)
1895 : :
1896 : : ! If istep is not equal to 1
1897 [ + - ]: 44 : else if(istep>=2)then
1898 : : !
1899 : : ! Compute the corresponding operator expectation value
1900 : : ! And put the residual vector minus the difference
1901 : : ! between vtrial and vtrial_old
1902 : : ! (this is actually the action of the operator !) in vect(*,2)
1903 : 44 : eigen_scf=0.0_dp
1904 [ + + ]: 102 : do isp=1,nspden
1905 [ + + ]: 524118 : do ifft=1,nfft
1906 : : eigen_scf=eigen_scf+&
1907 : 524074 : & work(ifft,isp,1) * vrespc(ifft,isp)
1908 : : end do
1909 : : end do
1910 : :
1911 [ + + ]: 102 : do isp=1,nspden
1912 [ + + ]: 524118 : do ifft=1,nfft
1913 : : vrespc(ifft,isp)=vrespc(ifft,isp)&
1914 : 524016 : & +vtrial(ifft,isp)-vtrial0(ifft,isp)
1915 : 524074 : work(ifft,isp,2)=vrespc(ifft,isp)
1916 : : end do
1917 : : end do
1918 : 44 : eigen_scf=eigen_scf/fix_resid
1919 : : write(message, '(a,es12.4,a)' ) &
1920 : 44 : & ' scfeig : Operator expectation value ',eigen_scf,' (extremal eigenvalue * diemix)'
1921 : 44 : call wrtout(std_out,message,'COLL')
1922 : 44 : call wrtout(ab_out,message,'COLL')
1923 : : !
1924 : : ! Compute residual of vect(*,2)
1925 : 44 : resid_new=zero
1926 [ + + ]: 102 : do isp=1,min(nspden,2)
1927 [ + + ]: 524118 : do ifft=1,nfft
1928 : 524074 : resid_new=resid_new+ work(ifft,isp,2) ** 2
1929 : : end do
1930 : : end do
1931 : : if (nspden==4) then
1932 : : do ifft=1,nfft
1933 : : resid_new=resid_new+two*(work(ifft,3,2)**2+work(ifft,4,2)**2)
1934 : : end do
1935 : : end if
1936 : 44 : factor=sqrt(fix_resid/resid_new)
1937 [ + - ]: 44 : if(eigen_scf<zero) then
1938 : 44 : factor=-factor ! the new vector MAY be oposite to the old one
1939 : : ! if(factor<-one) factor=-factor ! the new vector is not opposed to the old one
1940 : : end if
1941 : : write(message, '(a,es12.4)' ) &
1942 : 44 : & ' scfeig : Inverse of renormalization factor ',one/factor
1943 : 44 : call wrtout(std_out,message,'COLL')
1944 : 44 : call wrtout(ab_out,message,'COLL')
1945 : : write(message, '(a,es12.4)' ) &
1946 : 44 : & ' scfeig : Convergence criterion value (->0 at convergency) ',one/factor-eigen_scf-one
1947 : 44 : call wrtout(std_out,message,'COLL')
1948 : 44 : call wrtout(ab_out,message,'COLL')
1949 : :
1950 [ + + + + ]: 524118 : work(:,:,1)=work(:,:,2)*factor
1951 [ + + + + ]: 524118 : vtrial(:,:)=vtrial0(:,:)+work(:,:,1)
1952 : : ! End the different istep cases
1953 : : end if
1954 : :
1955 : : end subroutine scfeig
1956 : : !!***
1957 : :
1958 : : !!****f* m_abi_mixing/scfopt
1959 : : !!
1960 : : !! NAME
1961 : : !! scfopt
1962 : : !!
1963 : : !! FUNCTION
1964 : : !! Compute the next vtrial of the SCF cycle.
1965 : : !! Possible algorithms are : simple mixing, Anderson (order 1 or 2), Pulay
1966 : : !!
1967 : : !! INPUTS
1968 : : !! cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
1969 : : !! iscf= 2 => simple mixing
1970 : : !! = 3,4 => Anderson mixing
1971 : : !! = 7 => Pulay mixing
1972 : : !! istep= number of the step in the SCF cycle
1973 : : !! mpicomm=the mpi communicator used for the summation
1974 : : !! comm_atom=the mpi communicator over atoms ; PAW only (optional argument)
1975 : : !! mpi_summarize=set it to .true. if parallelisation is done over FFT
1976 : : !! nfft=(effective) number of FFT grid points (for this processor)
1977 : : !! npawmix=-PAW only- number of spherical part elements to be mixed
1978 : : !! nspden=number of spin-density components
1979 : : !! n_fftgr=third dimension of the array f_fftgr
1980 : : !! n_index=dimension for indices of potential/density (see ivrespc, i_vtrial...)
1981 : : !! opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential
1982 : : !! 1 vtrial (and also f_fftgr) actually contains the trial density
1983 : : !! pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part
1984 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
1985 : : !!
1986 : : !! OUTPUT
1987 : : !! (see side effects)
1988 : : !!
1989 : : !! SIDE EFFECTS
1990 : : !! vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave
1991 : : !! the input preconditioned residual potential
1992 : : !! at output, it is the new trial potential .
1993 : : !! f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid :
1994 : : !! The input vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(1)).
1995 : : !! The old vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(2)).
1996 : : !! The input preconditioned residual potential is in f_fftgr(:,:,i_vrespc(1))
1997 : : !! Two input old preconditioned residual potentials in f_fftgr(:,:,i_vrespc(2)) and f_fftgr(:,:,i_vrespc(3))
1998 : : !! Before output a permutation of i_vrespc(1), i_vrespc(2) and i_vrespc(3) occurs, without
1999 : : !! actually copying all the data (change of pointer).
2000 : : !! i_vrespc(n_index)=index of the preconditioned residual potentials (present and past) in the array f_fftgr
2001 : : !! i_vtrial(n_index) =indices of the potential (present and past) in the array f_fftgr
2002 : : !! ==== if usepaw==1
2003 : : !! f_paw(npawmix,n_fftgr*mffmem*usepaw)=different functions used for PAW
2004 : : !! (same as f_fftgr but for spherical part)
2005 : : !! vpaw(npawmix*usepaw)=at input, the aug. occupancies (rhoij) that gave
2006 : : !! the input preconditioned residual potential
2007 : : !! at output, it is the new aug. occupancies.
2008 : : !!
2009 : : !! SOURCE
2010 : :
2011 : 55039 : subroutine scfopt(cplex,f_fftgr,f_paw,iscf,istep,i_vrespc,i_vtrial,&
2012 : : & mpicomm,mpi_summarize,nfft,npawmix,nspden,n_fftgr,&
2013 : 55039 : & n_index,opt_denpot,pawoptmix,usepaw,vpaw,vresid,vtrial,errid,errmess, &
2014 : 55039 : & useextfpmd,f_extfpmd,nelect_extfpmd,&
2015 : 55039 : & use_rcpaw,nrcpawmix,f_rcpaw,rcpaw_occ,&
2016 : : & comm_atom) ! optional
2017 : :
2018 : : !Arguments ------------------------------------
2019 : : !scalars
2020 : : integer,intent(in) :: cplex,iscf,istep,n_fftgr,n_index,nfft,useextfpmd,use_rcpaw,nrcpawmix
2021 : : integer,intent(in) :: npawmix,nspden,opt_denpot,pawoptmix,usepaw,mpicomm
2022 : : integer, intent(in),optional :: comm_atom
2023 : : integer,intent(out) :: errid
2024 : : character(len = 500), intent(out) :: errmess
2025 : : logical, intent(in) :: mpi_summarize
2026 : : real(dp), intent(out) :: vresid
2027 : : real(dp),intent(inout) :: nelect_extfpmd
2028 : : real(dp),intent(inout) :: rcpaw_occ(nrcpawmix*use_rcpaw)
2029 : : !arrays
2030 : : integer,intent(inout) :: i_vrespc(n_index),i_vtrial(n_index)
2031 : : real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr)
2032 : : real(dp),intent(inout) :: f_paw(npawmix,n_fftgr*usepaw),vpaw(npawmix*usepaw)
2033 : : real(dp),intent(inout) :: vtrial(cplex*nfft,nspden)
2034 : : real(dp),intent(inout) :: f_extfpmd(useextfpmd*n_fftgr)
2035 : : real(dp),intent(inout) :: f_rcpaw(use_rcpaw*nrcpawmix,use_rcpaw*n_fftgr)
2036 : : !Local variables-------------------------------
2037 : : !scalars
2038 : : integer,parameter :: npulaymax=50
2039 : : integer :: i_vstore,ierr,ifft,ii,index,isp,jj,comm_atom_,niter,npulay,tmp
2040 : : real(dp),save :: prod_resid_old,resid_old,resid_old2
2041 : : real(dp) :: aa1,aa2,bb,cc1,cc2,current,det,lambda,lambda2,resid_best
2042 : : character(len=500) :: message
2043 : : !arrays
2044 : 55039 : integer,allocatable :: ipiv(:)
2045 : : real(dp),save :: amat(npulaymax+1,npulaymax+1)
2046 : : real(dp) :: mpibuff(2),prod_resid(1),prod_resid2(1),resid_new(1)
2047 : 55039 : real(dp),allocatable :: alpha(:),amatinv(:,:),amat_paw(:),rwork(:)
2048 : :
2049 : : ! *************************************************************************
2050 : :
2051 : : !FIXME Decreasing optimization with NVHPC in this routine
2052 : : ! because of significant numerical divergence (spotted in 24.9)
2053 : : !nvf$r opt=1
2054 : :
2055 : : !DEBUG
2056 : : !write(std_out,*)' scfopt : enter ; istep,iscf ',istep,iscf
2057 : : !ENDDEBUG
2058 : :
2059 : 55039 : errid = AB7_NO_ERROR
2060 : :
2061 [ + - ]: 55039 : comm_atom_=xmpi_comm_self; if(present(comm_atom)) comm_atom_=comm_atom
2062 : :
2063 : 55039 : i_vstore=i_vtrial(1)
2064 [ + + ]: 55039 : if (iscf==4) i_vstore=i_vtrial(2)
2065 [ + + ]: 55039 : if (iscf==7) then
2066 [ + - ]: 54528 : if (modulo(n_fftgr, 2) == 0 ) then
2067 : 54528 : npulay=(n_fftgr-2)/2
2068 : : else
2069 : 0 : npulay=(n_fftgr-1)/2
2070 : : end if
2071 : 54528 : i_vstore=i_vtrial(npulay)
2072 : : else
2073 : : npulay=0
2074 : : end if
2075 : :
2076 : : !Compute the new residual resid_new, from f_fftgr/f_paw(:,:,i_vrespc(1))
2077 : 55039 : call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr)
2078 [ + + + + ]: 55039 : if (usepaw==1.and.pawoptmix==1) then
2079 [ + + ]: 7680 : do index=1,npawmix
2080 : 7680 : resid_new(1)=resid_new(1)+f_paw(index,i_vrespc(1))**2
2081 : : end do
2082 : 160 : call xmpi_sum(resid_new(1),comm_atom_,ierr)
2083 : : end if
2084 : 55039 : vresid = resid_new(1)
2085 : :
2086 : : !_______________________________________________________________
2087 : : !Here use only the preconditioning, or initialize the other algorithms
2088 : :
2089 [ + + + + ]: 55039 : if (istep==1 .or. iscf==2) then
2090 : 8135 : write(message,'(2a)') ch10,' Simple mixing update:'
2091 : 8135 : call wrtout(std_out,message,'COLL')
2092 : :
2093 : 8135 : write(message,*)' residual square of the potential: ',resid_new(1)
2094 : 8135 : call wrtout(std_out,message,'COLL')
2095 : :
2096 : : ! Store information for later use
2097 [ + + ]: 8135 : if (iscf==3.or.iscf==4) resid_old=resid_new(1)
2098 [ + + ]: 8135 : if (iscf==7) then
2099 : 7849 : amat(:,:)=zero
2100 : 7849 : amat(1,1)=resid_new(1)
2101 : : end if
2102 : :
2103 : : ! Compute new vtrial (and new rhoij if PAW)
2104 [ + + + + : 182554740 : if (iscf/=2) f_fftgr(:,:,i_vstore)=vtrial(:,:)
+ + ]
2105 [ + + + + ]: 187770000 : vtrial(:,:)=vtrial(:,:)+f_fftgr(:,:,i_vrespc(1))
2106 [ + + ]: 8135 : if (usepaw==1) then
2107 [ + + + + ]: 385191 : if (iscf/=2) f_paw(:,i_vstore)=vpaw(:)
2108 [ + + ]: 386487 : vpaw(:)=vpaw(:)+f_paw(:,i_vrespc(1))
2109 : : end if
2110 : :
2111 [ + + ]: 8135 : if(useextfpmd==1) then
2112 [ + - ]: 8 : if (iscf/=2) f_extfpmd(i_vstore)=nelect_extfpmd
2113 : 8 : nelect_extfpmd=nelect_extfpmd+f_extfpmd(i_vrespc(1))
2114 : : endif
2115 : :
2116 [ + + ]: 8135 : if(use_rcpaw==1) then
2117 [ + - + + ]: 32 : if (iscf/=2) f_rcpaw(:,i_vstore)=rcpaw_occ(:)
2118 [ + + ]: 32 : rcpaw_occ(:)=rcpaw_occ(:)+f_rcpaw(:,i_vrespc(1))
2119 : : endif
2120 : :
2121 : : ! _______________________________________________________________
2122 : : ! Here Anderson algorithm using one previous iteration
2123 [ + + + + : 46904 : else if((istep==2 .or. iscf==3).and.iscf/=7)then
+ + ]
2124 : :
2125 : 202 : write(message,'(2a)') ch10,' Anderson update:'
2126 : 202 : call wrtout(std_out,message,'COLL')
2127 : :
2128 : 202 : write(message,*)' residual square of the potential: ',resid_new(1)
2129 : 202 : call wrtout(std_out,message,'COLL')
2130 : :
2131 : : ! Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2))
2132 : : call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,&
2133 : 202 : & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2134 [ + + - + ]: 202 : if (usepaw==1.and.pawoptmix==1) then
2135 [ # # ]: 0 : do index=1,npawmix
2136 : 0 : prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2))
2137 : : end do
2138 : 0 : call xmpi_sum(prod_resid(1),comm_atom_,ierr)
2139 : : end if
2140 : :
2141 : : ! Compute mixing factor
2142 : 202 : lambda=(resid_new(1)-prod_resid(1))/(resid_new(1)+resid_old-2*prod_resid(1))
2143 : 202 : write(message,*)' mixing of old trial potential: ',lambda
2144 : 202 : call wrtout(std_out,message,'COLL')
2145 : :
2146 : : ! Evaluate best residual square on the line
2147 : : resid_best=(1.0_dp-lambda)*(1.0_dp-lambda)*resid_new(1)&
2148 : : & +(1.0_dp-lambda)*lambda *2*prod_resid(1)&
2149 : 202 : & +lambda *lambda *resid_old
2150 : 202 : write(message,*)' predicted best residual square on the line: ',resid_best
2151 : 202 : call wrtout(std_out,message,'COLL')
2152 : :
2153 : : ! Store information for later use
2154 [ + + ]: 202 : if (iscf==4) then
2155 : 2 : prod_resid_old=prod_resid(1)
2156 : 2 : resid_old2=resid_old
2157 : : end if
2158 : 202 : resid_old=resid_new(1)
2159 : :
2160 : : ! Save latest trial potential and compute new trial potential
2161 [ + + ]: 490 : do isp=1,nspden
2162 [ + + ]: 2321642 : do ifft=1,cplex*nfft
2163 : 2321152 : current=vtrial(ifft,isp)
2164 : : vtrial(ifft,isp)=(one-lambda)*(current +f_fftgr(ifft,isp,i_vrespc(1)))&
2165 : 2321152 : & +lambda *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))
2166 : 2321440 : f_fftgr(ifft,isp,i_vstore)=current
2167 : : end do
2168 : : end do
2169 : :
2170 : : ! PAW: save latest rhoij and compute new rhoij
2171 [ + + ]: 43996 : do index=1,npawmix
2172 : 43794 : current=vpaw(index)
2173 : : vpaw(index)=(one-lambda)*(current +f_paw(index,i_vrespc(1)))&
2174 : 43794 : & +lambda *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))
2175 : 43996 : f_paw(index,i_vstore)=current
2176 : : end do
2177 : :
2178 [ - + ]: 202 : if(useextfpmd==1) then
2179 : 0 : current=nelect_extfpmd
2180 : 0 : nelect_extfpmd=(one-lambda)*f_extfpmd(i_vrespc(1))+lambda*f_extfpmd(i_vtrial(1))+f_extfpmd(i_vrespc(2))
2181 : 0 : f_extfpmd(i_vstore)=current
2182 : : endif
2183 : :
2184 [ + - ]: 202 : do index=1,nrcpawmix
2185 : 0 : current=rcpaw_occ(index)
2186 : : rcpaw_occ(index)=(one-lambda)*(current+f_rcpaw(index,i_vrespc(1)))&
2187 : 0 : & +lambda *(f_rcpaw(index,i_vtrial(1))+f_rcpaw(index,i_vrespc(2)))
2188 : 202 : f_rcpaw(index,i_vstore)=current
2189 : : end do
2190 : :
2191 : :
2192 : : ! _______________________________________________________________
2193 : : ! Here Anderson algorithm using two previous iterations
2194 [ + + ]: 46702 : else if(iscf==4.and.iscf/=7)then
2195 : :
2196 : 23 : write(message,'(2a)') ch10,' Anderson (order 2) update:'
2197 : 23 : call wrtout(std_out,message,'COLL')
2198 : :
2199 : 23 : write(message,*)' residual square of the potential: ',resid_new(1)
2200 : 23 : call wrtout(std_out,message,'COLL')
2201 : :
2202 : : ! Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2))
2203 : : call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,&
2204 : 23 : & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2205 [ + - - + ]: 23 : if (usepaw==1.and.pawoptmix==1) then
2206 [ # # ]: 0 : do index=1,npawmix
2207 : 0 : prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2))
2208 : : end do
2209 : : end if
2210 : :
2211 : : ! Compute prod_resid2 from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(3))
2212 : : call dotprodm_v(cplex,1,prod_resid2,i_vrespc(1),i_vrespc(3),mpicomm,mpi_summarize,1,1,&
2213 : 23 : & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2214 [ + - - + ]: 23 : if (usepaw==1.and.pawoptmix==1) then
2215 [ # # ]: 0 : do index=1,npawmix
2216 : 0 : prod_resid2(1)=prod_resid2(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(3))
2217 : : end do
2218 : : ! MPI reduction
2219 : 0 : mpibuff(1)=prod_resid(1);mpibuff(2)=prod_resid2(1)
2220 : 0 : call xmpi_sum(mpibuff,comm_atom_,ierr)
2221 : 0 : prod_resid(1)=mpibuff(1);prod_resid2(1)=mpibuff(2)
2222 : : end if
2223 : :
2224 : : ! Compute mixing factors
2225 : 23 : aa1=resid_new(1)+resid_old -two*prod_resid (1)
2226 : 23 : aa2=resid_new(1)+resid_old2-two*prod_resid2(1)
2227 : 23 : bb =resid_new(1)+prod_resid_old-prod_resid(1)-prod_resid2(1)
2228 : 23 : cc1=resid_new(1)-prod_resid (1)
2229 : 23 : cc2=resid_new(1)-prod_resid2(1)
2230 : 23 : det=aa1*aa2-bb*bb
2231 : 23 : lambda =(aa2*cc1-bb*cc2)/det
2232 : 23 : lambda2=(aa1*cc2-bb*cc1)/det
2233 : 23 : write(message,*)' mixing of old trial potentials: ',lambda,lambda2
2234 : 23 : call wrtout(std_out,message,'COLL')
2235 : :
2236 : : ! Store information for later use
2237 : 23 : prod_resid_old=prod_resid(1)
2238 : 23 : resid_old2=resid_old
2239 : 23 : resid_old=resid_new(1)
2240 : :
2241 : : ! Save latest trial potential and compute new trial potential
2242 [ + + ]: 46 : do isp=1,nspden
2243 [ + + ]: 88366 : do ifft=1,cplex*nfft
2244 : 88320 : current=vtrial(ifft,isp)
2245 : : vtrial(ifft,isp)=&
2246 : : & (one-lambda-lambda2)*(current +f_fftgr(ifft,isp,i_vrespc(1)))&
2247 : : & +lambda *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))&
2248 : 88320 : & +lambda2 *(f_fftgr(ifft,isp,i_vtrial(2))+f_fftgr(ifft,isp,i_vrespc(3)))
2249 : 88343 : f_fftgr(ifft,isp,i_vstore)=current
2250 : : end do
2251 : : end do
2252 : :
2253 : : ! PAW: save latest rhoij and compute new rhoij
2254 [ + + ]: 851 : do index=1,npawmix
2255 : 828 : current=vpaw(index)
2256 : : vpaw(index)=&
2257 : : & (one-lambda-lambda2)*(current +f_paw(index,i_vrespc(1)))&
2258 : : & +lambda *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))&
2259 : 828 : & +lambda2 *(f_paw(index,i_vtrial(2))+f_paw(index,i_vrespc(3)))
2260 : 851 : f_paw(index,i_vstore)=current
2261 : : end do
2262 : :
2263 [ - + ]: 23 : if(useextfpmd==1) then
2264 : 0 : current=nelect_extfpmd
2265 : : nelect_extfpmd=(one-lambda-lambda2)*(current+f_extfpmd(i_vrespc(1)))&
2266 : : & +lambda*(f_extfpmd(i_vtrial(1))+f_extfpmd(i_vrespc(2)))&
2267 : 0 : & +lambda2*(f_extfpmd(i_vtrial(2))+f_extfpmd(i_vrespc(3)))
2268 : 0 : f_extfpmd(i_vstore)=current
2269 : : endif
2270 : :
2271 [ + - ]: 23 : do index=1,nrcpawmix
2272 : 0 : current=rcpaw_occ(index)
2273 : : rcpaw_occ(index)=&
2274 : : & (one-lambda-lambda2)*(current +f_rcpaw(index,i_vrespc(1)))&
2275 : : & +lambda *(f_rcpaw(index,i_vtrial(1))+f_rcpaw(index,i_vrespc(2)))&
2276 : 0 : & +lambda2 *(f_rcpaw(index,i_vtrial(2))+f_rcpaw(index,i_vrespc(3)))
2277 : 23 : f_rcpaw(index,i_vstore)=current
2278 : : end do
2279 : :
2280 : :
2281 : : ! _______________________________________________________________
2282 : : ! Here Pulay algorithm
2283 [ + - ]: 46679 : else if(iscf==7)then
2284 : :
2285 : 46679 : niter=min(istep,npulay+1)
2286 : :
2287 : 46679 : write(message,'(2a,i2,a)') ch10,' Pulay update with ',niter-1,' previous iterations:'
2288 : 46679 : call wrtout(std_out,message,'COLL')
2289 : :
2290 [ - + ]: 46679 : if (npulay>npulaymax) then
2291 : 0 : errid = AB7_ERROR_MIXING_CONVERGENCE
2292 : 0 : write(errmess, '(4a)' ) ch10,&
2293 : 0 : & ' scfopt: ERROR - ',ch10,&
2294 : 0 : & ' Too many iterations required for Pulay algorithm (<50) !'
2295 : : return
2296 : : end if
2297 : :
2298 : : ! Compute "A" matrix
2299 [ + + ]: 46679 : if (istep>npulay+1) then
2300 [ + + ]: 85935 : do jj=1,niter-1
2301 [ + + ]: 613462 : do ii=1,niter-1
2302 : 602706 : amat(ii,jj)=amat(ii+1,jj+1)
2303 : : end do
2304 : : end do
2305 : : end if
2306 [ + + + + ]: 46679 : if (usepaw==1.and.pawoptmix==1) then
2307 [ - + + - : 429 : ABI_MALLOC(amat_paw,(niter))
+ - + - -
+ - + ]
2308 [ + + ]: 934 : amat_paw(:)=zero
2309 [ + + ]: 934 : do ii=1,niter
2310 [ + + ]: 37110 : do index=1,npawmix
2311 : 36967 : amat_paw(ii)=amat_paw(ii)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(1+niter-ii))
2312 : : end do
2313 : : end do
2314 : 143 : call xmpi_sum(amat_paw,comm_atom_,ierr)
2315 : : end if
2316 [ + + ]: 287634 : do ii=1,niter
2317 : : call dotprodm_v(cplex,1,amat(ii,niter),i_vrespc(1),i_vrespc(1+niter-ii),mpicomm,mpi_summarize,1,1,&
2318 : 240955 : & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2319 [ + + + + ]: 240955 : if (usepaw==1.and.pawoptmix==1) amat(ii,niter)=amat(ii,niter)+amat_paw(ii)
2320 [ + + ]: 287634 : if (ii<niter) amat(niter,ii)=amat(ii,niter)
2321 : : end do
2322 [ + + + + ]: 46679 : if (usepaw==1.and.pawoptmix==1)then
2323 [ - + ]: 143 : ABI_FREE(amat_paw)
2324 : : end if
2325 : :
2326 : : ! Invert "A" matrix
2327 [ - + + - : 186716 : ABI_MALLOC(amatinv,(niter,niter))
- + + - +
- + - - +
- + ]
2328 [ + + + + ]: 1769433 : amatinv(1:niter,1:niter)=amat(1:niter,1:niter)
2329 [ - + + - : 140037 : ABI_MALLOC(ipiv,(niter))
+ - + - -
+ - + ]
2330 [ - + + - : 140037 : ABI_MALLOC(rwork,(niter))
+ - + - -
+ - + ]
2331 : 46679 : call dgetrf(niter,niter,amatinv,niter,ipiv,ierr)
2332 : 46679 : call dgetri(niter,amatinv,niter,ipiv,rwork,niter,ierr)
2333 [ - + ]: 46679 : ABI_FREE(ipiv)
2334 [ - + ]: 46679 : ABI_FREE(rwork)
2335 : :
2336 : : ! Compute "alpha" factors
2337 [ - + + - : 140037 : ABI_MALLOC(alpha,(niter))
+ - + - -
+ - + ]
2338 : 46679 : det=zero
2339 [ + + ]: 287634 : alpha(:)=zero
2340 [ + + ]: 287634 : do ii=1,niter
2341 [ + + ]: 1769433 : do jj=1,niter
2342 : 1481799 : alpha(ii)=alpha(ii)+amatinv(jj,ii)
2343 : 1722754 : det=det+amatinv(jj,ii)
2344 : : end do
2345 : : end do
2346 [ + + ]: 287634 : alpha(:)=alpha(:)/det
2347 [ - + ]: 46679 : ABI_FREE(amatinv)
2348 : 46679 : write(message,'(a,5(1x,g10.3))')' mixing of old trial potential: alpha(m:m-4)=',(alpha(ii),ii=niter,max(1,niter-4),-1)
2349 : 46679 : call wrtout(std_out,message,'COLL')
2350 : :
2351 : : ! Save latest trial potential and compute new trial potential
2352 [ + + ]: 105632 : do isp=1,nspden
2353 [ + + ]: 1113371701 : do ifft=1,cplex*nfft
2354 : 1113266069 : current=vtrial(ifft,isp)
2355 : 1113266069 : vtrial(ifft,isp)=alpha(niter)*(current+f_fftgr(ifft,isp,i_vrespc(1)))
2356 [ + + ]: 6224245133 : do ii=niter-1,1,-1
2357 : : vtrial(ifft,isp)=vtrial(ifft,isp)+alpha(ii) &
2358 : 6224245133 : & *(f_fftgr(ifft,isp,i_vtrial(niter-ii))+f_fftgr(ifft,isp,i_vrespc(1+niter-ii)))
2359 : : end do
2360 : 1113325022 : f_fftgr(ifft,isp,i_vstore)=current
2361 : : end do
2362 : : end do
2363 : :
2364 : : ! PAW: save latest rhoij and compute new rhoij
2365 [ + + ]: 4054853 : do index=1,npawmix
2366 : 4008174 : current=vpaw(index)
2367 : 4008174 : vpaw(index)=alpha(niter)*(current+f_paw(index,i_vrespc(1)))
2368 [ + + ]: 25669130 : do ii=niter-1,1,-1
2369 : : vpaw(index)=vpaw(index)+alpha(ii) &
2370 : 25669130 : & *(f_paw(index,i_vtrial(niter-ii))+f_paw(index,i_vrespc(1+niter-ii)))
2371 : : end do
2372 : 4054853 : f_paw(index,i_vstore)=current
2373 : : end do
2374 : :
2375 [ + + ]: 46679 : if(useextfpmd==1) then
2376 : 75 : current=nelect_extfpmd
2377 : 75 : nelect_extfpmd=alpha(niter)*(current+f_extfpmd(i_vrespc(1)))
2378 [ + + ]: 456 : do ii=niter-1,1,-1
2379 : : nelect_extfpmd=nelect_extfpmd+alpha(ii)&
2380 : 456 : & *(f_extfpmd(i_vtrial(niter-ii))+f_extfpmd(i_vrespc(1+niter-ii)))
2381 : : enddo
2382 : 75 : f_extfpmd(i_vstore)=current
2383 : : endif
2384 : :
2385 [ + + ]: 46709 : do index=1,nrcpawmix
2386 : 30 : current=rcpaw_occ(index)
2387 : 30 : rcpaw_occ(index)=alpha(niter)*(current+f_rcpaw(index,i_vrespc(1)))
2388 [ + + ]: 84 : do ii=niter-1,1,-1
2389 : : rcpaw_occ(index)=rcpaw_occ(index)+alpha(ii) &
2390 : 84 : & *(f_rcpaw(index,i_vtrial(niter-ii))+f_rcpaw(index,i_vrespc(1+niter-ii)))
2391 : : end do
2392 : 46709 : f_rcpaw(index,i_vstore)=current
2393 : : end do
2394 : :
2395 : :
2396 [ - + ]: 46679 : ABI_FREE(alpha)
2397 : : ! _______________________________________________________________
2398 : : ! End of choice of optimization method
2399 : : end if
2400 : :
2401 : : !Permute potential indices
2402 [ + + ]: 55039 : if (iscf==3) then
2403 : 216 : tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
2404 [ + + ]: 54823 : else if (iscf==4) then
2405 : 27 : tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
2406 : 27 : tmp=i_vtrial(2) ; i_vtrial(2)=i_vtrial(1) ; i_vtrial(1)=tmp
2407 [ + + ]: 54796 : else if (iscf==7) then
2408 : 54528 : tmp=i_vtrial( npulay)
2409 [ + + ]: 381898 : do ii= npulay,2,-1
2410 : 381898 : i_vtrial(ii)=i_vtrial(ii-1)
2411 : : end do
2412 : 54528 : i_vtrial(1)=tmp
2413 : 54528 : tmp=i_vrespc(1+npulay)
2414 [ + + ]: 436426 : do ii=1+npulay,2,-1
2415 : 436426 : i_vrespc(ii)=i_vrespc(ii-1)
2416 : : end do
2417 : 54528 : i_vrespc(1)=tmp
2418 : : end if
2419 : :
2420 [ - + ]: 55039 : end subroutine scfopt
2421 : : !!***
2422 : :
2423 [ - + ]: 9 : subroutine scfopt_pulay_delta_sp(cplex,f_fftgr,f_respc_sp,f_trial_sp,f_trial_delta_i2, &
2424 [ - + - + : 18 : & f_trial_delta_scale,f_paw,istep,i_vrespc,i_vtrial,mpicomm,mpi_summarize,nfft, &
- + ]
2425 [ - + ]: 9 : & npawmix,nspden,n_fftgr,n_index,opt_denpot,pawoptmix,usepaw,vpaw,vresid,vtrial, &
2426 : 9 : & errid,errmess,useextfpmd,f_extfpmd,nelect_extfpmd,use_rcpaw,nrcpawmix,f_rcpaw, &
2427 : 9 : & rcpaw_occ,comm_atom)
2428 : :
2429 : : !Arguments ------------------------------------
2430 : : !scalars
2431 : : integer,intent(in) :: cplex,istep,n_fftgr,n_index,nfft,useextfpmd,use_rcpaw,nrcpawmix
2432 : : integer,intent(in) :: npawmix,nspden,opt_denpot,pawoptmix,usepaw,mpicomm,comm_atom
2433 : : integer,intent(out) :: errid
2434 : : character(len = 500), intent(out) :: errmess
2435 : : logical, intent(in) :: mpi_summarize
2436 : : real(dp), intent(out) :: vresid
2437 : : real(dp),intent(inout) :: nelect_extfpmd
2438 : : real(dp),intent(inout) :: rcpaw_occ(nrcpawmix*use_rcpaw)
2439 : : !arrays
2440 : : integer,intent(inout) :: i_vrespc(n_index),i_vtrial(n_index)
2441 : : real(dp),intent(inout) :: f_fftgr(:,:,:)
2442 : : real(sp),intent(inout) :: f_respc_sp(:,:,:),f_trial_sp(:,:)
2443 : : integer(c_int16_t),intent(inout) :: f_trial_delta_i2(:,:,:)
2444 : : real(dp),intent(inout) :: f_trial_delta_scale(:)
2445 : : real(dp),intent(inout) :: f_paw(npawmix,n_fftgr*usepaw),vpaw(npawmix*usepaw)
2446 : : real(dp),intent(inout) :: vtrial(cplex*nfft,nspden)
2447 : : real(dp),intent(inout) :: f_extfpmd(useextfpmd*n_fftgr)
2448 : : real(dp),intent(inout) :: f_rcpaw(use_rcpaw*nrcpawmix,use_rcpaw*n_fftgr)
2449 : : !Local variables-------------------------------
2450 : : !scalars
2451 : : integer,parameter :: npulaymax=50,pulay_delta_qmax=32767
2452 : : integer :: ierr,ifft,ii,index,isp,jj,kk,niter,npulay,tmp
2453 : : integer :: i_vstore,order,respc_comp,trial_comp,qint
2454 : : real(dp) :: alpha_sum,current,det,max_delta,old_newest,trial_value
2455 : : logical :: pulay_fallback
2456 : : character(len=500) :: message
2457 : : !arrays
2458 : 9 : integer,allocatable :: ipiv(:)
2459 : : real(dp),save :: amat(npulaymax+1,npulaymax+1)
2460 : : real(dp) :: resid_new(1)
2461 : 9 : real(dp),allocatable :: alpha(:),amatinv(:,:),amat_paw(:),rwork(:)
2462 : :
2463 : : ! *************************************************************************
2464 : :
2465 : 9 : errid = AB7_NO_ERROR
2466 [ + - ]: 9 : if (modulo(n_fftgr, 2) == 0 ) then
2467 : 9 : npulay=(n_fftgr-2)/2
2468 : : else
2469 : 0 : npulay=(n_fftgr-1)/2
2470 : : end if
2471 : 9 : i_vstore=i_vtrial(npulay)
2472 : 9 : niter=min(istep,npulay+1)
2473 : :
2474 : : call dotprodm_pulay_dp(cplex,1,resid_new,f_fftgr(:,:,1),f_fftgr(:,:,1), &
2475 [ - + - - : 9 : & mpicomm,mpi_summarize,nfft,nspden,opt_denpot)
- - - - -
+ - - - -
- - - + -
+ ]
2476 [ + - - + ]: 9 : if (usepaw==1.and.pawoptmix==1) then
2477 [ # # ]: 0 : do index=1,npawmix
2478 : 0 : resid_new(1)=resid_new(1)+f_paw(index,i_vrespc(1))**2
2479 : : end do
2480 : 0 : call xmpi_sum(resid_new(1),comm_atom,ierr)
2481 : : end if
2482 : 9 : vresid = resid_new(1)
2483 : :
2484 [ + + ]: 9 : if (istep==1) then
2485 : 1 : amat(:,:)=zero
2486 : 1 : amat(1,1)=resid_new(1)
2487 : : end if
2488 : :
2489 : 9 : write(message,'(2a,i2,a)') ch10,' Pulay delta-encoded update with ',niter-1,' previous iterations:'
2490 : 9 : call wrtout(std_out,message,'COLL')
2491 : :
2492 [ - + ]: 9 : if (npulay>npulaymax) then
2493 : 0 : errid = AB7_ERROR_MIXING_CONVERGENCE
2494 : 0 : write(errmess, '(4a)' ) ch10,&
2495 : 0 : & ' scfopt_pulay_delta_sp: ERROR - ',ch10,&
2496 : 0 : & ' Too many iterations required for Pulay algorithm (<50) !'
2497 : 0 : return
2498 : : end if
2499 : :
2500 [ + + ]: 9 : if (istep>npulay+1) then
2501 [ + + ]: 8 : do jj=1,niter-1
2502 [ + + ]: 57 : do ii=1,niter-1
2503 : 56 : amat(ii,jj)=amat(ii+1,jj+1)
2504 : : end do
2505 : : end do
2506 : : end if
2507 [ + - - + ]: 9 : if (usepaw==1.and.pawoptmix==1) then
2508 [ # # # # : 0 : ABI_MALLOC(amat_paw,(niter))
# # # # #
# # # ]
2509 [ # # ]: 0 : amat_paw(:)=zero
2510 [ # # ]: 0 : do ii=1,niter
2511 [ # # ]: 0 : if (ii==niter) then
2512 [ # # ]: 0 : do index=1,npawmix
2513 : 0 : amat_paw(ii)=amat_paw(ii)+f_paw(index,i_vrespc(1))**2
2514 : : end do
2515 : : else
2516 [ # # ]: 0 : do index=1,npawmix
2517 : 0 : amat_paw(ii)=amat_paw(ii)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(1+niter-ii))
2518 : : end do
2519 : : end if
2520 : : end do
2521 : 0 : call xmpi_sum(amat_paw,comm_atom,ierr)
2522 : : end if
2523 [ + + ]: 53 : do ii=1,niter
2524 [ + + ]: 44 : if (ii==niter) then
2525 : : call dotprodm_pulay_dp(cplex,1,amat(ii,niter),f_fftgr(:,:,1),f_fftgr(:,:,1), &
2526 [ - + - - : 9 : & mpicomm,mpi_summarize,nfft,nspden,opt_denpot)
- - - - -
+ - - - -
- - - + -
+ ]
2527 : : else
2528 [ + + ]: 35 : if (modulo(i_vrespc(1+niter-ii),2)==0) then
2529 : 28 : respc_comp=i_vrespc(1+niter-ii)/2
2530 : : else
2531 : : respc_comp=npulay+1
2532 : : end if
2533 [ + - - + ]: 35 : if (respc_comp<1 .or. respc_comp>npulay+1) then
2534 [ # # # # ]: 0 : if (usepaw==1.and.pawoptmix==1) then
2535 [ # # ]: 0 : ABI_FREE(amat_paw)
2536 : : end if
2537 : 0 : errid = AB7_ERROR_MIXING_ARG
2538 : 0 : write(errmess, '(4a)' ) ch10,&
2539 : 0 : & ' scfopt_pulay_delta_sp: ERROR - ',ch10,&
2540 : 0 : & ' Unexpected preconditioned-residual slot for compact Pulay history.'
2541 : 0 : return
2542 : : end if
2543 : : call dotprodm_pulay_sp(cplex,1,amat(ii,niter),f_fftgr(:,:,1), &
2544 [ - + - - : 35 : & f_respc_sp(:,:,respc_comp),mpicomm,mpi_summarize,nfft,nspden,opt_denpot)
- - - - -
+ - - - -
- - - + -
+ ]
2545 : : end if
2546 [ + - - + ]: 44 : if (usepaw==1.and.pawoptmix==1) amat(ii,niter)=amat(ii,niter)+amat_paw(ii)
2547 [ + + ]: 53 : if (ii<niter) amat(niter,ii)=amat(ii,niter)
2548 : : end do
2549 [ + - - + ]: 9 : if (usepaw==1.and.pawoptmix==1)then
2550 [ # # ]: 0 : ABI_FREE(amat_paw)
2551 : : end if
2552 : :
2553 [ - + + - : 27 : ABI_MALLOC(alpha,(niter))
+ - + - -
+ - + ]
2554 : 9 : pulay_fallback=.false.
2555 [ + + ]: 9 : if (niter==1) then
2556 : 1 : alpha(1)=one
2557 : : else
2558 [ - + + - : 32 : ABI_MALLOC(amatinv,(niter,niter))
- + + - +
- + - - +
- + ]
2559 [ + + + + ]: 318 : amatinv(1:niter,1:niter)=amat(1:niter,1:niter)
2560 [ - + + - : 24 : ABI_MALLOC(ipiv,(niter))
+ - + - -
+ - + ]
2561 [ - + + - : 16 : ABI_MALLOC(rwork,(niter))
- + - + ]
2562 : 8 : call dgetrf(niter,niter,amatinv,niter,ipiv,ierr)
2563 [ + - ]: 8 : if (ierr==0) call dgetri(niter,amatinv,niter,ipiv,rwork,niter,ierr)
2564 [ - + ]: 8 : ABI_FREE(ipiv)
2565 [ - + ]: 8 : ABI_FREE(rwork)
2566 [ - + ]: 8 : if (ierr/=0) then
2567 : : pulay_fallback=.true.
2568 : : else
2569 [ + + ]: 51 : alpha=zero
2570 : 8 : det=zero
2571 [ + + ]: 51 : do ii=1,niter
2572 [ + + ]: 318 : do jj=1,niter
2573 : 267 : alpha(ii)=alpha(ii)+amatinv(jj,ii)
2574 : 310 : det=det+amatinv(jj,ii)
2575 : : end do
2576 : : end do
2577 [ - + - + ]: 8 : if (det/=det .or. abs(det)<tiny(one)) then
2578 : : pulay_fallback=.true.
2579 : : else
2580 [ + + ]: 51 : alpha(:)=alpha(:)/det
2581 [ + + ]: 51 : alpha_sum=sum(abs(alpha(:)))
2582 [ - + - + ]: 8 : if (alpha_sum/=alpha_sum .or. alpha_sum>50.0_dp) pulay_fallback=.true.
2583 : : end if
2584 : : end if
2585 [ - + ]: 8 : ABI_FREE(amatinv)
2586 : : end if
2587 [ + - ]: 8 : if (pulay_fallback) then
2588 [ # # ]: 0 : ABI_FREE(alpha)
2589 : 0 : niter=1
2590 [ # # ]: 0 : ABI_MALLOC(alpha,(niter))
2591 : 0 : alpha(1)=one
2592 : : end if
2593 : 9 : write(message,'(a,5(1x,g10.3))')' mixing of old trial potential: alpha(m:m-4)=',(alpha(ii),ii=niter,max(1,niter-4),-1)
2594 : 9 : call wrtout(std_out,message,'COLL')
2595 : :
2596 [ + + ]: 9 : if (istep>1) then
2597 : : max_delta=zero
2598 [ + + ]: 16 : do isp=1,nspden
2599 [ + + ]: 65552 : do ifft=1,cplex*nfft
2600 : 65544 : max_delta=max(max_delta,abs(real(f_trial_sp(ifft,isp),dp)-vtrial(ifft,isp)))
2601 : : end do
2602 : : end do
2603 : 8 : trial_comp=(i_vtrial(1)+1)/2
2604 : 8 : f_trial_delta_scale(trial_comp)=max(max_delta/real(pulay_delta_qmax,dp),tiny(one))
2605 : : end if
2606 : :
2607 [ + + ]: 18 : do isp=1,nspden
2608 [ + + ]: 73746 : do ifft=1,cplex*nfft
2609 : 73728 : current=vtrial(ifft,isp)
2610 : 73728 : vtrial(ifft,isp)=alpha(niter)*(current+f_fftgr(ifft,isp,1))
2611 [ + + ]: 360448 : do ii=niter-1,1,-1
2612 : 286720 : order=niter-ii
2613 : 286720 : old_newest=real(f_trial_sp(ifft,isp),dp)
2614 : 286720 : trial_value=old_newest
2615 [ + + ]: 917504 : do kk=2,order
2616 : 630784 : trial_comp=(i_vtrial(kk)+1)/2
2617 : 917504 : trial_value=trial_value+real(f_trial_delta_i2(ifft,isp,trial_comp),dp)*f_trial_delta_scale(trial_comp)
2618 : : end do
2619 [ + + ]: 286720 : if (modulo(i_vrespc(1+niter-ii),2)==0) then
2620 : 229376 : respc_comp=i_vrespc(1+niter-ii)/2
2621 : : else
2622 : : respc_comp=npulay+1
2623 : : end if
2624 : : vtrial(ifft,isp)=vtrial(ifft,isp)+alpha(ii) * &
2625 : 360448 : & (trial_value+real(f_respc_sp(ifft,isp,respc_comp),dp))
2626 : : end do
2627 [ + + ]: 73728 : if (istep>1) then
2628 : 65536 : trial_comp=(i_vtrial(1)+1)/2
2629 : 65536 : qint=nint((real(f_trial_sp(ifft,isp),dp)-current)/f_trial_delta_scale(trial_comp))
2630 : 65536 : qint=max(-pulay_delta_qmax,min(pulay_delta_qmax,qint))
2631 : 65536 : f_trial_delta_i2(ifft,isp,trial_comp)=int(qint,c_int16_t)
2632 : : end if
2633 : 73728 : f_trial_sp(ifft,isp)=real(current,sp)
2634 [ + + ]: 73728 : if (modulo(i_vrespc(1),2)==0) then
2635 : 65536 : respc_comp=i_vrespc(1)/2
2636 : : else
2637 : : respc_comp=npulay+1
2638 : : end if
2639 : 73737 : f_respc_sp(ifft,isp,respc_comp)=real(f_fftgr(ifft,isp,1),sp)
2640 : : end do
2641 : : end do
2642 : :
2643 [ + + ]: 657 : do index=1,npawmix
2644 : 648 : current=vpaw(index)
2645 : 648 : vpaw(index)=alpha(niter)*(current+f_paw(index,i_vrespc(1)))
2646 [ + + ]: 3168 : do ii=niter-1,1,-1
2647 : : vpaw(index)=vpaw(index)+alpha(ii) &
2648 : 3168 : & *(f_paw(index,i_vtrial(niter-ii))+f_paw(index,i_vrespc(1+niter-ii)))
2649 : : end do
2650 : 657 : f_paw(index,i_vstore)=current
2651 : : end do
2652 : :
2653 [ - + ]: 9 : if(useextfpmd==1) then
2654 : 0 : current=nelect_extfpmd
2655 : 0 : nelect_extfpmd=alpha(niter)*(current+f_extfpmd(i_vrespc(1)))
2656 [ # # ]: 0 : do ii=niter-1,1,-1
2657 : : nelect_extfpmd=nelect_extfpmd+alpha(ii)&
2658 : 0 : & *(f_extfpmd(i_vtrial(niter-ii))+f_extfpmd(i_vrespc(1+niter-ii)))
2659 : : enddo
2660 : 0 : f_extfpmd(i_vstore)=current
2661 : : endif
2662 : :
2663 [ - + ]: 9 : do index=1,nrcpawmix
2664 : 0 : current=rcpaw_occ(index)
2665 : 0 : rcpaw_occ(index)=alpha(niter)*(current+f_rcpaw(index,i_vrespc(1)))
2666 [ # # ]: 0 : do ii=niter-1,1,-1
2667 : : rcpaw_occ(index)=rcpaw_occ(index)+alpha(ii) &
2668 : 0 : & *(f_rcpaw(index,i_vtrial(niter-ii))+f_rcpaw(index,i_vrespc(1+niter-ii)))
2669 : : end do
2670 : 9 : f_rcpaw(index,i_vstore)=current
2671 : : end do
2672 : :
2673 [ - + ]: 9 : ABI_FREE(alpha)
2674 : :
2675 : 9 : tmp=i_vtrial(npulay)
2676 [ + + ]: 63 : do ii=npulay,2,-1
2677 : 63 : i_vtrial(ii)=i_vtrial(ii-1)
2678 : : end do
2679 : 9 : i_vtrial(1)=tmp
2680 : 9 : tmp=i_vrespc(1+npulay)
2681 [ + + ]: 72 : do ii=1+npulay,2,-1
2682 : 72 : i_vrespc(ii)=i_vrespc(ii-1)
2683 : : end do
2684 : 9 : i_vrespc(1)=tmp
2685 : :
2686 [ - + - - ]: 9 : end subroutine scfopt_pulay_delta_sp
2687 : : !!***
2688 : :
2689 : 18 : subroutine dotprodm_pulay_dp(cplex,cpldot,dot,potarr1,potarr2,mpicomm,mpi_summarize,nfft,nspden,opt_storage)
2690 : :
2691 : : !Arguments ------------------------------------
2692 : : integer,intent(in) :: cpldot,cplex,nfft,nspden,opt_storage,mpicomm
2693 : : logical, intent(in) :: mpi_summarize
2694 : : real(dp),intent(in) :: potarr1(cplex*nfft,nspden)
2695 : : real(dp),intent(in) :: potarr2(cplex*nfft,nspden)
2696 : : real(dp),intent(out) :: dot(cpldot)
2697 : :
2698 : : !Local variables-------------------------------
2699 : : integer :: ierr,ifft,ispden
2700 : : real(dp) :: ai,ar
2701 : : ! *************************************************************************
2702 : :
2703 : : DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
2704 : : DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot")
2705 : : DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
2706 : :
2707 [ + - + - ]: 18 : if(cplex==1 .or. cpldot==1)then
2708 : 18 : ar=zero
2709 [ + + ]: 36 : do ispden=1,min(nspden,2)
2710 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
2711 [ + + ]: 147492 : do ifft=1,cplex*nfft
2712 : 147474 : ar=ar + potarr1(ifft,ispden)*potarr2(ifft,ispden)
2713 : : end do
2714 : : end do
2715 : 18 : dot(1)=ar
2716 [ - + ]: 18 : if (nspden==4) then
2717 : : ar=zero
2718 [ # # ]: 0 : do ispden=3,4
2719 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
2720 [ # # ]: 0 : do ifft=1,cplex*nfft
2721 : 0 : ar=ar + potarr1(ifft,ispden)*potarr2(ifft,ispden)
2722 : : end do
2723 : : end do
2724 [ # # ]: 0 : if (opt_storage==0) then
2725 [ # # ]: 0 : if (cplex==1) then
2726 : 0 : dot(1)=dot(1)+two*ar
2727 : : else
2728 : 0 : dot(1)=dot(1)+ar
2729 : : end if
2730 : : else
2731 : 0 : dot(1)=half*(dot(1)+ar)
2732 : : end if
2733 : : end if
2734 : : else
2735 : 0 : ar=zero ; ai=zero
2736 [ # # ]: 0 : do ispden=1,min(nspden,2)
2737 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
2738 [ # # ]: 0 : do ifft=1,nfft
2739 : : ar=ar + potarr1(2*ifft-1,ispden)*potarr2(2*ifft-1,ispden) &
2740 : 0 : & + potarr1(2*ifft ,ispden)*potarr2(2*ifft ,ispden)
2741 : : ai=ai + potarr1(2*ifft-1,ispden)*potarr2(2*ifft ,ispden) &
2742 : 0 : & - potarr1(2*ifft ,ispden)*potarr2(2*ifft-1,ispden)
2743 : : end do
2744 : : end do
2745 : 0 : dot(1)=ar ; dot(2)=ai
2746 [ # # ]: 0 : if (nspden==4) then
2747 : : ar=zero
2748 [ # # ]: 0 : do ispden=3,4
2749 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
2750 [ # # ]: 0 : do ifft=1,nfft
2751 : : ar=ar + potarr1(2*ifft-1,ispden)*potarr2(2*ifft-1,ispden) &
2752 : 0 : & + potarr1(2*ifft ,ispden)*potarr2(2*ifft ,ispden)
2753 : : ai=ai + potarr1(2*ifft-1,ispden)*potarr2(2*ifft ,ispden) &
2754 : 0 : & - potarr1(2*ifft ,ispden)*potarr2(2*ifft-1,ispden)
2755 : : end do
2756 : : end do
2757 [ # # ]: 0 : if (opt_storage==0) then
2758 : 0 : dot(1)=dot(1)+ar
2759 : 0 : dot(2)=dot(2)+ai
2760 : : else
2761 : 0 : dot(1)=half*(dot(1)+ar)
2762 : 0 : dot(2)=half*(dot(2)+ai)
2763 : : end if
2764 : : end if
2765 : : end if
2766 : :
2767 [ - + ]: 18 : if (mpi_summarize) call xmpi_sum(dot,mpicomm,ierr)
2768 [ - + - - ]: 18 : if(cpldot==2 .and. cplex==1)dot(2)=zero
2769 : :
2770 : 18 : end subroutine dotprodm_pulay_dp
2771 : : !!***
2772 : :
2773 : : !!****f* ABINIT/dotprodm_pulay_sp
2774 : : !! NAME
2775 : : !! dotprodm_pulay_sp
2776 : : !!
2777 : : !! FUNCTION
2778 : : !! Double-precision dot product between a double-precision current vector
2779 : : !! and a single-precision historical vector.
2780 : : !!
2781 : : !! SOURCE
2782 : :
2783 : 35 : subroutine dotprodm_pulay_sp(cplex,cpldot,dot,potarr1,potarr2,mpicomm,mpi_summarize,nfft,nspden,opt_storage)
2784 : :
2785 : : !Arguments ------------------------------------
2786 : : integer,intent(in) :: cpldot,cplex,nfft,nspden,opt_storage,mpicomm
2787 : : logical, intent(in) :: mpi_summarize
2788 : : real(dp),intent(in) :: potarr1(cplex*nfft,nspden)
2789 : : real(sp),intent(in) :: potarr2(cplex*nfft,nspden)
2790 : : real(dp),intent(out) :: dot(cpldot)
2791 : :
2792 : : !Local variables-------------------------------
2793 : : integer :: ierr,ifft,ispden
2794 : : real(dp) :: ai,ar
2795 : : ! *************************************************************************
2796 : :
2797 : : DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
2798 : : DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot")
2799 : : DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
2800 : :
2801 [ + - + - ]: 35 : if(cplex==1 .or. cpldot==1)then
2802 : 35 : ar=zero
2803 [ + + ]: 70 : do ispden=1,min(nspden,2)
2804 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
2805 [ + + ]: 286790 : do ifft=1,cplex*nfft
2806 : 286755 : ar=ar + potarr1(ifft,ispden)*real(potarr2(ifft,ispden),dp)
2807 : : end do
2808 : : end do
2809 : 35 : dot(1)=ar
2810 [ - + ]: 35 : if (nspden==4) then
2811 : : ar=zero
2812 [ # # ]: 0 : do ispden=3,4
2813 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
2814 [ # # ]: 0 : do ifft=1,cplex*nfft
2815 : 0 : ar=ar + potarr1(ifft,ispden)*real(potarr2(ifft,ispden),dp)
2816 : : end do
2817 : : end do
2818 [ # # ]: 0 : if (opt_storage==0) then
2819 [ # # ]: 0 : if (cplex==1) then
2820 : 0 : dot(1)=dot(1)+two*ar
2821 : : else
2822 : 0 : dot(1)=dot(1)+ar
2823 : : end if
2824 : : else
2825 : 0 : dot(1)=half*(dot(1)+ar)
2826 : : end if
2827 : : end if
2828 : : else
2829 : 0 : ar=zero ; ai=zero
2830 [ # # ]: 0 : do ispden=1,min(nspden,2)
2831 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
2832 [ # # ]: 0 : do ifft=1,nfft
2833 : : ar=ar + potarr1(2*ifft-1,ispden)*real(potarr2(2*ifft-1,ispden),dp) &
2834 : 0 : & + potarr1(2*ifft ,ispden)*real(potarr2(2*ifft ,ispden),dp)
2835 : : ai=ai + potarr1(2*ifft-1,ispden)*real(potarr2(2*ifft ,ispden),dp) &
2836 : 0 : & - potarr1(2*ifft ,ispden)*real(potarr2(2*ifft-1,ispden),dp)
2837 : : end do
2838 : : end do
2839 : 0 : dot(1)=ar ; dot(2)=ai
2840 [ # # ]: 0 : if (nspden==4) then
2841 : : ar=zero
2842 [ # # ]: 0 : do ispden=3,4
2843 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
2844 [ # # ]: 0 : do ifft=1,nfft
2845 : : ar=ar + potarr1(2*ifft-1,ispden)*real(potarr2(2*ifft-1,ispden),dp) &
2846 : 0 : & + potarr1(2*ifft ,ispden)*real(potarr2(2*ifft ,ispden),dp)
2847 : : ai=ai + potarr1(2*ifft-1,ispden)*real(potarr2(2*ifft ,ispden),dp) &
2848 : 0 : & - potarr1(2*ifft ,ispden)*real(potarr2(2*ifft-1,ispden),dp)
2849 : : end do
2850 : : end do
2851 [ # # ]: 0 : if (opt_storage==0) then
2852 : 0 : dot(1)=dot(1)+ar
2853 : 0 : dot(2)=dot(2)+ai
2854 : : else
2855 : 0 : dot(1)=half*(dot(1)+ar)
2856 : 0 : dot(2)=half*(dot(2)+ai)
2857 : : end if
2858 : : end if
2859 : : end if
2860 : :
2861 [ - + ]: 35 : if (mpi_summarize) call xmpi_sum(dot,mpicomm,ierr)
2862 [ - + - - ]: 35 : if(cpldot==2 .and. cplex==1)dot(2)=zero
2863 : :
2864 : 35 : end subroutine dotprodm_pulay_sp
2865 : : !!***
2866 : :
2867 : : !!****f* ABINIT/findminscf
2868 : : !! NAME
2869 : : !! findminscf
2870 : : !!
2871 : : !! FUNCTION
2872 : : !! Compute the minimum of a function whose value
2873 : : !! and derivative are known at two points, using different algorithms.
2874 : : !! Also deduce different quantities at this predicted
2875 : : !! point, and at the two other points
2876 : : !!
2877 : : !! INPUTS
2878 : : !! choice=1,uses a linear interpolation of the derivatives
2879 : : !! =2,uses a quadratic interpolation based on the
2880 : : !! values of the function, and the second derivative at mid-point
2881 : : !! etotal_1=first value of the function
2882 : : !! etotal_2=second value of the function
2883 : : !! dedv_1=first value of the derivative
2884 : : !! dedv_2=second value of the derivative
2885 : : !! lambda_1=first value of the argument
2886 : : !! lambda_2=second value of the argument
2887 : : !!
2888 : : !! OUTPUT
2889 : : !! dedv_predict=predicted value of the derivative (usually zero,
2890 : : !! except if choice=4, if it happens that a minimum cannot be located,
2891 : : !! and a trial step is taken)
2892 : : !! d2edv2_predict=predicted value of the second derivative (not if choice=4)
2893 : : !! d2edv2_1=first value of the second derivative (not if choice=4)
2894 : : !! d2edv2_2=second value of the second derivative (not if choice=4)
2895 : : !! etotal_predict=predicted value of the function
2896 : : !! lambda_predict=predicted value of the argument
2897 : : !! status= 0 if everything went normally ;
2898 : : !! 1 if negative second derivative
2899 : : !! 2 if some other problem
2900 : : !!
2901 : : !! SOURCE
2902 : :
2903 : 91 : subroutine findminscf(choice,dedv_1,dedv_2,dedv_predict,&
2904 : : & d2edv2_1,d2edv2_2,d2edv2_predict,&
2905 : : & etotal_1,etotal_2,etotal_predict,&
2906 : : & lambda_1,lambda_2,lambda_predict,errid,errmess)
2907 : :
2908 : : !Arguments ------------------------------------
2909 : : !scalars
2910 : : integer,intent(in) :: choice
2911 : : integer,intent(out) :: errid
2912 : : character(len=500), intent(out) :: errmess
2913 : : real(dp),intent(in) :: dedv_1,dedv_2,etotal_1,etotal_2,lambda_1,lambda_2
2914 : : real(dp),intent(out) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_predict
2915 : : real(dp),intent(out) :: etotal_predict,lambda_predict
2916 : :
2917 : : !Local variables-------------------------------
2918 : : !scalars
2919 : : real(dp) :: cc,d2edv2_mid,d_lambda,dedv_2bis
2920 : : real(dp) :: dedv_mid2,etotal_2bis
2921 : : character(len=500) :: message
2922 : :
2923 : : ! *************************************************************************
2924 : :
2925 : : !DEBUG
2926 : : !write(std_out,*)' findmin : enter'
2927 : : !write(std_out,*)' choice,lambda_1,lambda_2=',choice,lambda_1,lambda_2
2928 : : !ENDDEBUG
2929 : :
2930 : 91 : errid = AB7_NO_ERROR
2931 : 91 : d_lambda=lambda_1-lambda_2
2932 : :
2933 [ + + ]: 91 : if(choice==1) then
2934 : :
2935 : : ! Use the derivative information to predict lambda
2936 : 84 : d2edv2_mid=(dedv_1-dedv_2)/d_lambda
2937 : 84 : lambda_predict=lambda_2-dedv_2/d2edv2_mid
2938 : 84 : dedv_predict=dedv_2+(lambda_predict-lambda_2)*d2edv2_mid
2939 : 84 : d2edv2_1=d2edv2_mid
2940 : 84 : d2edv2_2=d2edv2_mid
2941 : 84 : d2edv2_predict=d2edv2_mid
2942 : : ! also use the first energy to predict new energy
2943 : : etotal_predict=etotal_1+dedv_1*(lambda_predict-lambda_1)&
2944 : 84 : & +0.5_dp*d2edv2_1*(lambda_predict-lambda_1)**2
2945 : : etotal_2bis=etotal_1+dedv_1*(lambda_2-lambda_1)&
2946 : 84 : & +0.5_dp*d2edv2_1*(lambda_2-lambda_1)**2
2947 : :
2948 [ + + ]: 84 : if(d2edv2_mid<0.0_dp)then
2949 : 4 : errid = AB7_ERROR_MIXING_INTERNAL
2950 : 4 : write(errmess,'(a,es18.10,a)')'The second derivative is negative, equal to ',d2edv2_mid,'.'
2951 : 4 : ABI_WARNING(errmess)
2952 : : end if
2953 : :
2954 [ + - ]: 7 : else if(choice==2) then
2955 : :
2956 : : ! Use energies and first derivative information
2957 : : ! etotal = aa + bb * lambda + cc * lambda**2
2958 : 7 : dedv_mid2=(etotal_1-etotal_2)/d_lambda
2959 : 7 : cc=(dedv_1-dedv_mid2)/d_lambda
2960 : 7 : lambda_predict=lambda_1-0.5_dp*dedv_1/cc
2961 : 7 : d2edv2_1=2*cc
2962 : 7 : d2edv2_2=d2edv2_1
2963 : 7 : d2edv2_predict=d2edv2_1
2964 [ - + ]: 7 : if(d2edv2_predict<0.0_dp)then
2965 : 0 : errid = AB7_ERROR_MIXING_INTERNAL
2966 : : write(errmess, '(a,es18.10,a,a,a)' )&
2967 : 0 : & 'The second derivative is negative, equal to',d2edv2_predict,'.',ch10,&
2968 : 0 : & '=> Pivoting '
2969 : 0 : ABI_WARNING(errmess)
2970 [ # # ]: 0 : if(etotal_2 < etotal_1)then
2971 : 0 : lambda_predict=lambda_2-0.5_dp*(lambda_1-lambda_2)
2972 : : else
2973 : 0 : lambda_predict=lambda_1-0.5_dp*(lambda_2-lambda_1)
2974 : : end if
2975 : : end if
2976 : 7 : dedv_predict=dedv_1+(lambda_predict-lambda_1)*d2edv2_1
2977 : 7 : dedv_2bis=dedv_1+(lambda_2-lambda_1)*d2edv2_1
2978 : : etotal_predict=etotal_1+dedv_1*(lambda_predict-lambda_1)&
2979 : 7 : & +0.5_dp*d2edv2_1*(lambda_predict-lambda_1)**2
2980 : :
2981 : : end if
2982 : :
2983 : : write(message, '(a,es12.4,a,es18.10)' ) &
2984 : 91 : & ' findmin : lambda_predict ',lambda_predict,' etotal_predict ',etotal_predict
2985 : 91 : call wrtout(std_out,message,'COLL')
2986 : :
2987 : 91 : end subroutine findminscf
2988 : : !!***
2989 : :
2990 : : !!****f* ABINIT/dotprodm_v
2991 : : !! NAME
2992 : : !! dotprodm_v
2993 : : !!
2994 : : !! FUNCTION
2995 : : !! For two sets of potentials,
2996 : : !! compute dot product of each pair of two potentials (integral over FFT grid), to obtain
2997 : : !! a series of square residual-like quantity (so the sum of product of values
2998 : : !! is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume).
2999 : : !! Take into account the spin components of the potentials (nspden),
3000 : : !! and sum over them.
3001 : : !! Need the index of the first pair of potentials to be treated, in each array
3002 : : !! of potentials, and the number of potentials to be treated.
3003 : : !! Might be used to compute just one square of norm, in
3004 : : !! a big array, such as to avoid copying a potential from a big array
3005 : : !! to a temporary place.
3006 : : !!
3007 : : !! INPUTS
3008 : : !! cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
3009 : : !! cpldot=if 1, the dot array is real, if 2, the dot array is complex
3010 : : !! index1=index of the first potential to be treated in the potarr1 array
3011 : : !! index2=index of the first potential to be treated in the potarr2 array
3012 : : !! mpicomm=the mpi communicator used for the summation
3013 : : !! mpi_summarize=set it to .true. if parallelisation is done over FFT
3014 : : !! mult1=number of potentials to be treated in the first set
3015 : : !! mult2=number of potentials to be treated in the second set
3016 : : !! nfft= (effective) number of FFT grid points (for this processor)
3017 : : !! npot1= third dimension of the potarr1 array
3018 : : !! npot2= third dimension of the potarr2 array
3019 : : !! nspden=number of spin-density components
3020 : : !! opt_storage: 0, if potentials are stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn]
3021 : : !! 1, if potentials are stored as V, B_x, B_y, Bz (B=magn. field)
3022 : : !! potarr1(cplex*nfft,nspden,npot)=first array of real space potentials on FFT grid
3023 : : !! (if cplex=2 and cpldot=2, potarr1 is the array that will be complex conjugated)
3024 : : !! potarr2(cplex*nfft,nspden,npot)=second array of real space potentials on FFT grid
3025 : : !!
3026 : : !! OUTPUT
3027 : : !! dot(cpldot,mult1,mult2)= series of values of the dot product
3028 : : !!
3029 : : !! SIDE EFFECTS
3030 : : !!
3031 : : !! NOTES
3032 : : !! Concerning storage when nspden=4:
3033 : : !! cplex=1:
3034 : : !! opt_storage=0: V are stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
3035 : : !! opt_storage=1: V are stored as : V, B_x, B_y, B_z (real)
3036 : : !! cplex=2:
3037 : : !! opt_storage=0: V are stored as : V^11, V^22, V^12, i.V^21 (complex)
3038 : : !! opt_storage=1: V are stored as : V, B_x, B_y, B_z (complex)
3039 : : !!
3040 : : !! SOURCE
3041 : :
3042 : 241203 : subroutine dotprodm_v(cplex,cpldot,dot,index1,index2,mpicomm,mpi_summarize,&
3043 : 241203 : & mult1,mult2,nfft,npot1,npot2,nspden,opt_storage,potarr1,potarr2)
3044 : :
3045 : : !Arguments ------------------------------------
3046 : : !scalars
3047 : : integer,intent(in) :: cpldot,cplex,index1,index2,mult1,mult2,nfft,npot1,npot2
3048 : : integer,intent(in) :: nspden,opt_storage,mpicomm
3049 : : logical, intent(in) :: mpi_summarize
3050 : : !arrays
3051 : : real(dp),intent(in) :: potarr1(cplex*nfft,nspden,npot1)
3052 : : real(dp),intent(in) :: potarr2(cplex*nfft,nspden,npot2)
3053 : : real(dp),intent(out) :: dot(cpldot,mult1,mult2)
3054 : :
3055 : : !Local variables-------------------------------
3056 : : !scalars
3057 : : integer :: i1,i2,ierr,ifft,ispden
3058 : : real(dp) :: ai,ar
3059 : : !arrays
3060 : : real(dp) :: tsec(2)
3061 : :
3062 : : ! *************************************************************************
3063 : :
3064 : : !Real or complex inputs are coded
3065 : : DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
3066 : :
3067 : : !Real or complex outputs are coded
3068 : : DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot")
3069 : : DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
3070 : : DBG_CHECK( npot1-index1-mult1 >= -1,"npot1-index1-mult1")
3071 : : DBG_CHECK( npot2-index2-mult2 >= -1,"npot2-index2-mult2")
3072 : :
3073 [ + + - + ]: 241203 : if(cplex==1 .or. cpldot==1)then
3074 : :
3075 [ + + ]: 482406 : do i1=1,mult1
3076 [ + + ]: 723609 : do i2=1,mult2
3077 : 241203 : ar=zero
3078 [ + + ]: 529205 : do ispden=1,min(nspden,2)
3079 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
3080 [ + + ]: 5892228922 : do ifft=1,cplex*nfft
3081 : 5891987719 : ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1)
3082 : : end do
3083 : : end do
3084 : 241203 : dot(1,i1,i2)=ar
3085 [ + + ]: 482406 : if (nspden==4) then
3086 : : ar=zero
3087 [ + + ]: 49482 : do ispden=3,4
3088 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
3089 [ + + ]: 335092690 : do ifft=1,cplex*nfft
3090 : 335076196 : ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1)
3091 : : end do
3092 : : end do
3093 [ + + ]: 16494 : if (opt_storage==0) then
3094 [ + + ]: 9756 : if (cplex==1) then
3095 : 9208 : dot(1,i1,i2)=dot(1,i1,i2)+two*ar
3096 : : else
3097 : 548 : dot(1,i1,i2)=dot(1,i1,i2)+ar
3098 : : end if
3099 : : else
3100 : 6738 : dot(1,i1,i2)=half*(dot(1,i1,i2)+ar)
3101 : : end if
3102 : : end if
3103 : : end do
3104 : : end do
3105 : :
3106 : : else ! if (cplex==2 .and. cpldot==2)
3107 : :
3108 [ # # ]: 0 : do i1=1,mult1
3109 [ # # ]: 0 : do i2=1,mult2
3110 : 0 : ar=zero ; ai=zero
3111 [ # # ]: 0 : do ispden=1,min(nspden,2)
3112 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
3113 [ # # ]: 0 : do ifft=1,nfft
3114 : : ar=ar + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) &
3115 : 0 : & + potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1)
3116 : : ai=ai + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1) &
3117 : 0 : & - potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1)
3118 : : end do
3119 : : end do
3120 : 0 : dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai
3121 [ # # ]: 0 : if (nspden==4) then
3122 : : ar=zero
3123 [ # # ]: 0 : do ispden=3,4
3124 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
3125 [ # # ]: 0 : do ifft=1,nfft
3126 : : ar=ar + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) &
3127 : 0 : & + potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1)
3128 : : ai=ai + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft ,ispden,index2+i2-1) &
3129 : 0 : & - potarr1(2*ifft ,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1)
3130 : : end do
3131 : : end do
3132 [ # # ]: 0 : if (opt_storage==0) then
3133 : 0 : dot(1,i1,i2)=dot(1,i1,i2)+ar
3134 : 0 : dot(2,i1,i2)=dot(2,i1,i2)+ai
3135 : : else
3136 : 0 : dot(1,i1,i2)=half*(dot(1,i1,i2)+ar)
3137 : 0 : dot(2,i1,i2)=half*(dot(2,i1,i2)+ai)
3138 : : end if
3139 : : end if
3140 : : end do
3141 : : end do
3142 : : end if
3143 : :
3144 : : !XG030513 : MPIWF reduction (addition) on dot is needed here
3145 [ + + ]: 241203 : if (mpi_summarize) then
3146 : 32428 : call timab(48,1,tsec)
3147 : 32428 : call xmpi_sum(dot,mpicomm ,ierr)
3148 : 32428 : call timab(48,2,tsec)
3149 : : end if
3150 : :
3151 [ - + - - : 241203 : if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero
- - - - ]
3152 : :
3153 : 241203 : end subroutine dotprodm_v
3154 : : !!***
3155 : :
3156 : : !!****f* ABINIT/dotprodm_vn
3157 : : !! NAME
3158 : : !! dotprodm_vn
3159 : : !!
3160 : : !! FUNCTION
3161 : : !! For a set of densities and a set of potentials,
3162 : : !! compute the dot product (integral over FFT grid) of each pair, to obtain
3163 : : !! a series of energy-like quantity (so the usual dotproduct is divided
3164 : : !! by the number of FFT points, and multiplied by the primitive cell volume).
3165 : : !! Take into account the spin components of the density and potentials (nspden),
3166 : : !! and sum correctly over them. Note that the storage of densities and
3167 : : !! potentials is different : for potential, one stores the matrix components,
3168 : : !! while for the density, one stores the trace, and then, either the
3169 : : !! spin-polarisation (if nspden=2), or the magnetization vector (if nspden=4).
3170 : : !! Need the index of the first density/potential pair to be treated, in each array,
3171 : : !! and the number of pairs to be treated.
3172 : : !! Might be used to compute just one dot product, in
3173 : : !! a big array, such as to avoid copying the density and potential from a big array
3174 : : !! to a temporary place.
3175 : : !!
3176 : : !! INPUTS
3177 : : !! cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
3178 : : !! cpldot=if 1, the dot array is real, if 2, the dot array is complex (not coded yet for nspden=4)
3179 : : !! denarr(cplex*nfft,nspden,nden)=real space density on FFT grid
3180 : : !! id=index of the first density to be treated in the denarr array
3181 : : !! ip=index of the first potential to be treated in the potarr array
3182 : : !! mpicomm=the mpi communicator used for the summation
3183 : : !! mpi_summarize=set it to .true. if parallelisation is done over FFT
3184 : : !! multd=number of densities to be treated
3185 : : !! multp=number of potentials to be treated
3186 : : !! nden=third dimension of the denarr array
3187 : : !! nfft= (effective) number of FFT grid points (for this processor)
3188 : : !! nfftot= total number of FFT grid points
3189 : : !! npot=third dimension of the potarr array
3190 : : !! nspden=number of spin-density components
3191 : : !! potarr(cplex*nfft,nspden,npot)=real space potential on FFT grid
3192 : : !! (will be complex conjugated if cplex=2 and cpldot=2)
3193 : : !! ucvol=unit cell volume (Bohr**3)
3194 : : !!
3195 : : !! OUTPUT
3196 : : !! dot(cpldot,multp,multd)= series of values of the dot product potential/density
3197 : : !!
3198 : : !! SIDE EFFECTS
3199 : : !!
3200 : : !! NOTES
3201 : : !! Concerning storage when nspden=4:
3202 : : !! cplex=1:
3203 : : !! V are stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
3204 : : !! N are stored as : n, m_x, m_y, m_z (real)
3205 : : !! cplex=2:
3206 : : !! V are stored as : V^11, V^22, V^12, i.V^21 (complex)
3207 : : !! N are stored as : n, m_x, m_y, mZ (complex)
3208 : : !!
3209 : : !! SOURCE
3210 : :
3211 : 377 : subroutine dotprodm_vn(cplex,cpldot,denarr,dot,id,ip,mpicomm, mpi_summarize,multd,multp,&
3212 : 377 : & nden,nfft,nfftot,npot,nspden,potarr,ucvol)
3213 : :
3214 : : !Arguments ------------------------------------
3215 : : !scalars
3216 : : integer,intent(in) :: cpldot,cplex,id,ip,multd,multp,nden,nfft,nfftot,npot
3217 : : integer,intent(in) :: nspden,mpicomm
3218 : : logical, intent(in) :: mpi_summarize
3219 : : real(dp),intent(in) :: ucvol
3220 : : !arrays
3221 : : real(dp),intent(in) :: denarr(cplex*nfft,nspden,nden)
3222 : : real(dp),intent(in) :: potarr(cplex*nfft,nspden,npot)
3223 : : real(dp),intent(out) :: dot(cpldot,multp,multd)
3224 : :
3225 : : !Local variables-------------------------------
3226 : : !scalars
3227 : : integer :: i1,i2,ierr,ir,jr
3228 : : real(dp) :: ai,ar,dim11,dim12,dim21,dim22,dim_dn,dim_up,dre11,dre12,dre21
3229 : : real(dp) :: dre22,dre_dn,dre_up,factor,pim11,pim12,pim21,pim22,pim_dn,pim_up
3230 : : real(dp) :: pre11,pre12,pre21,pre22,pre_dn,pre_up
3231 : : !arrays
3232 : : real(dp) :: tsec(2)
3233 : :
3234 : : ! *************************************************************************
3235 : :
3236 : : !Real or complex inputs are coded
3237 : : DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
3238 : :
3239 : : !Real or complex outputs are coded
3240 : : DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot")
3241 : : DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
3242 : :
3243 : : DBG_CHECK(id >= 1,'Wrong id')
3244 : : DBG_CHECK(ip >= 1,'Wrong id')
3245 : :
3246 : : DBG_CHECK(multd >= 1,"wrong multd")
3247 : : DBG_CHECK(multp >= 1,"wrong multp")
3248 : :
3249 : : DBG_CHECK(nden-id-multd >=-1,'nden-id-multd')
3250 : : DBG_CHECK(npot-ip-multp >=-1,'npot-ip-multp')
3251 : :
3252 [ + + ]: 377 : if(nspden==1)then
3253 : :
3254 [ - + - - ]: 257 : if(cpldot==1 .or. cplex==1 )then
3255 : :
3256 [ + + ]: 514 : do i2=1,multd
3257 [ + + ]: 771 : do i1=1,multp
3258 : : ar=zero
3259 : : !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar)
3260 [ + + ]: 994791 : do ir=1,cplex*nfft
3261 : 994791 : ar=ar + potarr(ir,1,ip+i1-1)*denarr(ir,1,id+i2-1)
3262 : : end do
3263 : 514 : dot(1,i1,i2)=ar
3264 : : end do ! i1
3265 : : end do ! i2
3266 : :
3267 : : else ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot
3268 : :
3269 [ # # ]: 0 : do i2=1,multd
3270 [ # # ]: 0 : do i1=1,multp
3271 : : ar=zero ; ai=zero
3272 : : !$OMP PARALLEL DO PRIVATE(ir,jr) SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai)
3273 [ # # ]: 0 : do ir=1,nfft
3274 : 0 : jr=2*ir
3275 : : ar=ar + potarr(jr-1,1,ip+i1-1)*denarr(jr-1,1,id+i2-1) &
3276 : 0 : & + potarr(jr ,1,ip+i1-1)*denarr(jr ,1,id+i2-1)
3277 : : ai=ai + potarr(jr-1,1,ip+i1-1)*denarr(jr ,1,id+i2-1) &
3278 : 0 : & - potarr(jr ,1,ip+i1-1)*denarr(jr-1,1,id+i2-1)
3279 : : end do
3280 : 0 : dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai
3281 : : end do ! i1
3282 : : end do ! i2
3283 : :
3284 : : end if
3285 : :
3286 [ + - ]: 120 : else if(nspden==2)then
3287 : :
3288 [ - + - - ]: 120 : if(cpldot==1 .or. cplex==1 )then
3289 : :
3290 [ + + ]: 240 : do i2=1,multd
3291 [ + + ]: 360 : do i1=1,multp
3292 : : ar=zero
3293 : : !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar)
3294 [ + + ]: 1624920 : do ir=1,cplex*nfft
3295 : : ar=ar + potarr(ir,1,ip+i1-1)* denarr(ir,2,id+i2-1) & ! This is the spin up contribution
3296 : 1624920 : & + potarr(ir,2,ip+i1-1)*(denarr(ir,1,id+i2-1)-denarr(ir,2,id+i2-1)) ! This is the spin down contribution
3297 : : end do
3298 : 240 : dot(1,i1,i2)=ar
3299 : : end do ! i1
3300 : : end do ! i2
3301 : :
3302 : : else ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot
3303 : :
3304 [ # # ]: 0 : do i2=1,multd
3305 [ # # ]: 0 : do i1=1,multp
3306 : : ar=zero ; ai=zero
3307 : : !$OMP PARALLEL DO PRIVATE(ir,jr,dre_up,dim_up,dre_dn,dim_dn,pre_up,pim_up,pre_dn,pim_dn) &
3308 : : !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai)
3309 [ # # ]: 0 : do ir=1,nfft
3310 : 0 : jr=2*ir
3311 : :
3312 : 0 : dre_up=denarr(jr-1,2,id+i2-1)
3313 : 0 : dim_up=denarr(jr ,2,id+i2-1)
3314 : 0 : dre_dn=denarr(jr-1,1,id+i2-1)-dre_up
3315 : 0 : dim_dn=denarr(jr ,1,id+i2-1)-dim_up
3316 : :
3317 : 0 : pre_up=potarr(jr-1,1,ip+i1-1)
3318 : 0 : pim_up=potarr(jr ,1,ip+i1-1)
3319 : 0 : pre_dn=potarr(jr-1,2,ip+i1-1)
3320 : 0 : pim_dn=potarr(jr ,2,ip+i1-1)
3321 : :
3322 : : ar=ar + pre_up * dre_up &
3323 : : & + pim_up * dim_up &
3324 : : & + pre_dn * dre_dn &
3325 : 0 : & + pim_dn * dim_dn
3326 : : ai=ai + pre_up * dim_up &
3327 : : & - pim_up * dre_up &
3328 : : & + pre_dn * dim_dn &
3329 : 0 : & - pim_dn * dre_dn
3330 : :
3331 : : end do
3332 : 0 : dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai
3333 : : end do ! i1
3334 : : end do ! i2
3335 : :
3336 : : end if
3337 : :
3338 [ # # ]: 0 : else if(nspden==4)then
3339 : : ! \rho{\alpha,\beta} V^{\alpha,\beta} =
3340 : : ! rho*(V^{11}+V^{22})/2$
3341 : : ! + m_x Re(V^{12})- m_y Im{V^{12}}+ m_z(V^{11}-V^{22})/2
3342 [ # # ]: 0 : if (cplex==1) then
3343 [ # # ]: 0 : do i2=1,multd
3344 [ # # ]: 0 : do i1=1,multp
3345 : : ar=zero
3346 : : !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar)
3347 [ # # ]: 0 : do ir=1,cplex*nfft
3348 : : ar=ar+(potarr(ir,1,ip+i1-1)+potarr(ir,2,ip+i1-1))*half*denarr(ir,1,id+i2-1)& ! This is the density contrib
3349 : : & + potarr(ir,3,ip+i1-1) *denarr(ir,2,id+i2-1)& ! This is the m_x contrib
3350 : : & - potarr(ir,4,ip+i1-1) *denarr(ir,3,id+i2-1)& ! This is the m_y contrib
3351 : 0 : & +(potarr(ir,1,ip+i1-1)-potarr(ir,2,ip+i1-1))*half*denarr(ir,4,id+i2-1) ! This is the m_z contrib
3352 : : end do
3353 : 0 : dot(1,i1,i2)=ar
3354 : : end do ! i1
3355 : : end do ! i2
3356 : : else ! cplex=2
3357 : : ! Note concerning storage when cplex=2:
3358 : : ! V are stored as : v^11, v^22, V^12, i.V^21 (each are complex)
3359 : : ! N are stored as : n, m_x, m_y, mZ (each are complex)
3360 [ # # ]: 0 : if (cpldot==1) then
3361 [ # # ]: 0 : do i2=1,multd
3362 [ # # ]: 0 : do i1=1,multp
3363 : 0 : ar=zero ; ai=zero
3364 : : !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre22,pim22,pre12,pim12) &
3365 : : !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar)
3366 [ # # ]: 0 : do ir=1,nfft
3367 : 0 : jr=2*ir
3368 : 0 : dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2))
3369 : 0 : dim11=half*(denarr(jr ,1,id+i2)+denarr(jr-1,4,id+i2))
3370 : 0 : dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2))
3371 : 0 : dim22=half*(denarr(jr ,1,id+i2)-denarr(jr-1,4,id+i2))
3372 : 0 : dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr ,3,id+i2))
3373 : 0 : dim12=half*(denarr(jr ,2,id+i2)-denarr(jr-1,3,id+i2))
3374 : 0 : dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr ,3,id+i2))
3375 : 0 : dim21=half*(denarr(jr ,2,id+i2)+denarr(jr-1,3,id+i2))
3376 : 0 : pre11= potarr(jr-1,1,ip+i1)
3377 : 0 : pim11= potarr(jr ,1,ip+i1)
3378 : 0 : pre22= potarr(jr-1,2,ip+i1)
3379 : 0 : pim22= potarr(jr ,2,ip+i1)
3380 : 0 : pre12= potarr(jr-1,3,ip+i1)
3381 : 0 : pim12= potarr(jr ,3,ip+i1)
3382 : 0 : pre21= potarr(jr ,4,ip+i1)
3383 : 0 : pim21=-potarr(jr-1,4,ip+i1)
3384 : : ar=ar + pre11 * dre11 &
3385 : : & + pim11 * dim11 &
3386 : : & + pre22 * dre22 &
3387 : : & + pim22 * dim22 &
3388 : : & + pre12 * dre12 &
3389 : : & + pim12 * dim12 &
3390 : : & + pre21 * dre21 &
3391 : 0 : & + pim21 * dim21
3392 : : end do
3393 : 0 : dot(1,i1,i2)=ar
3394 : : end do ! i1
3395 : : end do ! i2
3396 : : else !cpldot=2
3397 [ # # ]: 0 : do i2=1,multd
3398 [ # # ]: 0 : do i1=1,multp
3399 : : ar=zero ; ai=zero
3400 : : !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre12,pim12,pre22,pim22) &
3401 : : !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai)
3402 [ # # ]: 0 : do ir=1,nfft
3403 : 0 : jr=2*ir
3404 : 0 : dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2))
3405 : 0 : dim11=half*(denarr(jr ,1,id+i2)+denarr(jr-1,4,id+i2))
3406 : 0 : dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2))
3407 : 0 : dim22=half*(denarr(jr ,1,id+i2)-denarr(jr-1,4,id+i2))
3408 : 0 : dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr ,3,id+i2))
3409 : 0 : dim12=half*(denarr(jr ,2,id+i2)-denarr(jr-1,3,id+i2))
3410 : 0 : dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr ,3,id+i2))
3411 : 0 : dim21=half*(denarr(jr ,2,id+i2)+denarr(jr-1,3,id+i2))
3412 : 0 : pre11= potarr(jr-1,1,ip+i1)
3413 : 0 : pim11= potarr(jr ,1,ip+i1)
3414 : 0 : pre22= potarr(jr-1,2,ip+i1)
3415 : 0 : pim22= potarr(jr ,2,ip+i1)
3416 : 0 : pre12= potarr(jr-1,3,ip+i1)
3417 : 0 : pim12= potarr(jr ,3,ip+i1)
3418 : 0 : pre21= potarr(jr ,4,ip+i1)
3419 : 0 : pim21=-potarr(jr-1,4,ip+i1)
3420 : : ar=ar + pre11 * dre11 &
3421 : : & + pim11 * dim11 &
3422 : : & + pre22 * dre22 &
3423 : : & + pim22 * dim22 &
3424 : : & + pre12 * dre12 &
3425 : : & + pim12 * dim12 &
3426 : : & + pre21 * dre21 &
3427 : 0 : & + pim21 * dim21
3428 : : ai=ai + pre11 * dim11 &
3429 : : & - pim11 * dre11 &
3430 : : & + pre22 * dim22 &
3431 : : & - pim22 * dre22 &
3432 : : & + pre12 * dim12 &
3433 : : & - pim12 * dre12 &
3434 : : & + pre21 * dim21 &
3435 : 0 : & - pim21 * dre21
3436 : : end do
3437 : 0 : dot(1,i1,i2)=ar
3438 : 0 : dot(2,i1,i2)=ai
3439 : : end do ! i1
3440 : : end do ! i2
3441 : : end if ! cpldot
3442 : : end if ! cplex
3443 : : end if ! nspden
3444 : :
3445 : 377 : factor=ucvol/dble(nfftot)
3446 [ + + + + : 1508 : dot(:,:,:)=factor*dot(:,:,:)
+ + ]
3447 : :
3448 : : !XG030513 : MPIWF reduction (addition) on dot is needed here
3449 [ - + ]: 377 : if (mpi_summarize) then
3450 : 0 : call timab(48,1,tsec)
3451 : 0 : call xmpi_sum(dot,mpicomm ,ierr)
3452 : 0 : call timab(48,2,tsec)
3453 : : end if
3454 : :
3455 [ - + - - : 377 : if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero
- - - - ]
3456 : :
3457 : 377 : end subroutine dotprodm_vn
3458 : : !!***
3459 : :
3460 : : !!****f* ABINIT/sqnormm_v
3461 : : !! NAME
3462 : : !! sqnormm_v
3463 : : !!
3464 : : !! FUNCTION
3465 : : !! For a series of potentials,
3466 : : !! compute square of the norm (integral over FFT grid), to obtain
3467 : : !! a square residual-like quantity (so the sum of product of values
3468 : : !! is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume).
3469 : : !! Take into account the spin components of the density and potentials (nspden), and sum over them.
3470 : : !! Need the index of the first potential to be treated, in the provided array
3471 : : !! of potentials, and the number of potentials to be treated.
3472 : : !! Might be used to compute just one square of norm, in a big array, such as to avoid
3473 : : !! copying a potential from a big array to a temporary place.
3474 : : !!
3475 : : !! INPUTS
3476 : : !! cplex=if 1, real space function on FFT grid is REAL, if 2, COMPLEX
3477 : : !! index=index of the first potential to be treated
3478 : : !! mpicomm=the mpi communicator used for the summation
3479 : : !! mpi_summarize=set it to .true. if parallelisation is done over FFT
3480 : : !! mult=number of potentials to be treated
3481 : : !! nfft= (effective) number of FFT grid points (for this processor)
3482 : : !! npot= third dimension of the potarr array
3483 : : !! nspden=number of spin-density components
3484 : : !! opt_storage: 0, if potential is stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn]
3485 : : !! 1, if potential is stored as V, B_x, B_y, Bz (B=magn. field)
3486 : : !! potarr(cplex*nfft,nspden,npot)=array of real space potentials on FFT grid
3487 : : !!
3488 : : !! OUTPUT
3489 : : !! norm2(mult)= value of the square of the norm of the different potentials
3490 : : !!
3491 : : !! SOURCE
3492 : :
3493 : 55225 : subroutine sqnormm_v(cplex,index,mpicomm, mpi_summarize,mult,nfft,norm2,npot,nspden,opt_storage,potarr)
3494 : :
3495 : : !Arguments ------------------------------------
3496 : : !scalars
3497 : : integer,intent(in) :: cplex,index,mult,nfft,npot,nspden,opt_storage,mpicomm
3498 : : logical, intent(in) :: mpi_summarize
3499 : : !arrays
3500 : : real(dp),intent(in) :: potarr(cplex*nfft,nspden,npot)
3501 : : real(dp),intent(out) :: norm2(mult)
3502 : :
3503 : : !Local variables-------------------------------
3504 : : !scalars
3505 : : integer :: ierr,ifft,ii,ispden
3506 : : real(dp) :: ar
3507 : : !arrays
3508 : : real(dp) :: tsec(2)
3509 : :
3510 : : ! *************************************************************************
3511 : :
3512 : : !Real or complex inputs are coded
3513 : : DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
3514 : : DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
3515 : :
3516 : : DBG_CHECK(index>=1,"wrong index")
3517 : : DBG_CHECK(mult>=1,"wrong mult")
3518 : : DBG_CHECK(npot>=1,"wrong npot")
3519 : :
3520 : : DBG_CHECK(npot-index-mult>=-1,'npot-index-mult')
3521 : :
3522 [ + + ]: 110450 : do ii=1,mult
3523 : 55225 : ar=zero
3524 [ + + ]: 119000 : do ispden=1,min(nspden,2)
3525 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar)
3526 [ + + ]: 1252106039 : do ifft=1,cplex*nfft
3527 : 1252050814 : ar=ar + potarr(ifft,ispden,index+ii-1)**2
3528 : : end do
3529 : : end do
3530 : 55225 : norm2(ii)=ar
3531 [ + + ]: 110450 : if (nspden==4) then
3532 : : ar=zero
3533 [ + + ]: 7704 : do ispden=3,4
3534 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar)
3535 [ + + ]: 53516420 : do ifft=1,cplex*nfft
3536 : 53513852 : ar=ar + potarr(ifft,ispden,index+ii-1)**2
3537 : : end do
3538 : : end do
3539 [ + + ]: 2568 : if (opt_storage==0) then
3540 [ + + ]: 1577 : if (cplex==1) then
3541 : 1472 : norm2(ii)=norm2(ii)+two*ar
3542 : : else
3543 : 105 : norm2(ii)=norm2(ii)+ar
3544 : : end if
3545 : : else
3546 : 991 : norm2(ii)=half*(norm2(ii)+ar)
3547 : : end if
3548 : : end if
3549 : : end do
3550 : :
3551 : : !XG030513 : MPIWF reduction (addition) on norm2 is needed here
3552 [ + + ]: 55225 : if (mpi_summarize) then
3553 : 6876 : call timab(48,1,tsec)
3554 : 6876 : call xmpi_sum(norm2,mpicomm ,ierr)
3555 : 6876 : call timab(48,2,tsec)
3556 : : end if
3557 : :
3558 : 55225 : end subroutine sqnormm_v
3559 : : !!***
3560 : :
3561 : : !!****f* ABINIT/aprxdr
3562 : : !! NAME
3563 : : !! aprxdr
3564 : : !!
3565 : : !! FUNCTION
3566 : : !! Compute the approximative derivatives of the energy at different
3567 : : !! points along the line search, thanks to a finite-difference formula.
3568 : : !! This formula is the projection along the line search of the
3569 : : !! Eq.(11) in PRB54, 4383 (1996) [[cite:Gonze1996]].
3570 : : !!
3571 : : !! INPUTS
3572 : : !! cplex: if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
3573 : : !! choice= if==3, compute dedv_new, dedv_old, and dedv_mix,
3574 : : !! if/=3, compute only dedv_new and dedv_old.
3575 : : !! i_vresid and i_rhor, see the next lines.
3576 : : !! f_fftgr(nfft,nspden,n_fftgr)=different functions defined on the fft grid :
3577 : : !! The last residual potential is in f_fftgr(:,:,i_vresid(1)).
3578 : : !! The old residual potential is in f_fftgr(:,:,i_vresid(2)).
3579 : : !! The previous old residual potential is in f_fftgr(:,:,i_vresid(3)).
3580 : : !! (needed only when choice==3)
3581 : : !! The old density is in f_fftgr(:,:,i_rhor2).
3582 : : !! f_atm(3,natom,n_fftgr)=different functions defined for each atom :
3583 : : !! The last HF force is in f_atm(:,:,i_vresid(1)).
3584 : : !! The old HF force is in f_fftgr(:,:,i_vresid(2)).
3585 : : !! The previous old HF force is in f_fftgr(:,:,i_vresid(3)).
3586 : : !! (needed only when choice==3)
3587 : : !! The old atomic positions are in f_atm(:,:,i_rhor2)
3588 : : !! moved_atm_inside: if==1, the atoms are allowed to move.
3589 : : !! mpicomm=the mpi communicator used for the summation
3590 : : !! mpi_summarize=set it to .true. if parallelisation is done over FFT
3591 : : !! natom=number of atoms in unit cell
3592 : : !! nfft=(effective) number of FFT grid points (for this processor)
3593 : : !! nfftot=total number of FFT grid points
3594 : : !! nspden=number of spin-density components
3595 : : !! rhor(nfft,nspden)=actual density
3596 : : !! xred(3,natom)=reduced atomic coordinates
3597 : : !!
3598 : : !! OUTPUT
3599 : : !! dedv_mix=approximate derivative from previous old residual
3600 : : !! dedv_new=approximate derivative from new residual
3601 : : !! dedv_old=approximate derivative from old residual (output only when choice==3)
3602 : : !!
3603 : : !! NOTES
3604 : : !! Should be OpenMP parallelized
3605 : : !!
3606 : : !! SOURCE
3607 : :
3608 : 156 : subroutine aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,&
3609 : 156 : & f_atm,f_fftgr,i_rhor2,i_vresid,moved_atm_inside,&
3610 : 156 : & mpicomm,mpi_summarize,natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred)
3611 : :
3612 : : !Arguments ------------------------------------
3613 : : !scalars
3614 : : integer,intent(in) :: choice,cplex,i_rhor2,moved_atm_inside,n_fftgr,natom,nfft
3615 : : integer,intent(in) :: mpicomm,nfftot,nspden
3616 : : logical, intent(in) :: mpi_summarize
3617 : : real(dp),intent(in) :: ucvol
3618 : : real(dp),intent(out) :: dedv_mix,dedv_new,dedv_old
3619 : : !arrays
3620 : : integer,intent(in) :: i_vresid(3)
3621 : : real(dp),intent(in) :: f_atm(3,natom,n_fftgr)
3622 : : real(dp),intent(in) :: f_fftgr(cplex*nfft,nspden,n_fftgr)
3623 : : real(dp),intent(in) :: rhor(cplex*nfft,nspden),xred(3,natom)
3624 : :
3625 : : !Local variables-------------------------------
3626 : : !scalars
3627 : : integer :: iatom,idir
3628 : : !arrays
3629 : : real(dp) :: dedv_temp(1)
3630 : 156 : real(dp),allocatable :: ddens(:,:,:)
3631 : :
3632 : : ! *************************************************************************
3633 : :
3634 [ - + + - : 624 : ABI_MALLOC(ddens,(cplex*nfft,nspden,1))
- + + - +
- + - - +
- + ]
3635 : :
3636 : : !Compute approximative derivative of the energy
3637 : : !with respect to change of potential
3638 : :
3639 [ + + + + ]: 1757324 : ddens(:,:,1)=rhor(:,:)-f_fftgr(:,:,i_rhor2)
3640 : :
3641 : : !call dotprod_vn(cplex,1,ddens,dedv_old,nfft,nfftot,nspden,1,vresid,ucvol)
3642 : : !Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(2))
3643 : : call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(2),mpicomm,mpi_summarize,1,1,1,&
3644 : 156 : & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol)
3645 : 156 : dedv_old = dedv_temp(1)
3646 : :
3647 : : !Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(1))
3648 : : call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(1),mpicomm,mpi_summarize,1,1,1,&
3649 : 156 : & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol)
3650 : 156 : dedv_new= dedv_temp(1)
3651 : :
3652 [ + + ]: 156 : if(choice==3)then
3653 : : ! Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(3))
3654 : : call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(3),mpicomm,mpi_summarize,1,1,1,&
3655 : 65 : & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol)
3656 : 65 : dedv_mix = dedv_temp(1)
3657 : : end if
3658 : :
3659 [ - + ]: 156 : ABI_FREE(ddens)
3660 : :
3661 : : !-------------------------------------------------------
3662 : :
3663 : : !Now, take care of eventual atomic displacements
3664 : :
3665 [ + + ]: 156 : if(moved_atm_inside==1)then
3666 [ + + ]: 68 : do idir=1,3
3667 [ + + ]: 170 : do iatom=1,natom
3668 : : dedv_new=dedv_new+&
3669 : 102 : & f_atm(idir,iatom,i_vresid(1))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2))
3670 : : dedv_old=dedv_old+&
3671 : 102 : & f_atm(idir,iatom,i_vresid(2))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2))
3672 [ + + ]: 102 : if(choice==3) dedv_mix=dedv_mix+&
3673 : 87 : & f_atm(idir,iatom,i_vresid(3))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2))
3674 : : end do
3675 : : end do
3676 : : end if
3677 : :
3678 : 156 : end subroutine aprxdr
3679 : : !!***
3680 : :
3681 : 0 : end module m_abi_mixing
3682 : : !!***
|
