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 8080 : subroutine init_(mix)
111 :
112 : !Arguments ------------------------------------
113 : !scalars
114 : type(abi_mixing_object), intent(out) :: mix
115 : ! *************************************************************************
116 :
117 : ! Default values.
118 8080 : mix%iscf = ABI_MIXING_NONE
119 8080 : mix%mffmem = 1
120 8080 : mix%n_index = 0
121 8080 : mix%n_fftgr = 0
122 8080 : mix%n_pulayit = 7
123 8080 : mix%pulayhist_storage = ABI_MIXING_PULAY_STORAGE_FULL
124 8080 : mix%n_pawmix = 0
125 8080 : mix%n_atom = 0
126 8080 : mix%space = 0
127 8080 : mix%useprec = .true.
128 8080 : mix%useextfpmd = 0
129 8080 : mix%use_rcpaw = 0
130 8080 : mix%n_rcpawmix = 0
131 :
132 8080 : call nullify_(mix)
133 :
134 8080 : 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 16160 : subroutine nullify_(mix)
147 :
148 : !Arguments ------------------------------------
149 : !scalars
150 : type(abi_mixing_object), intent(inout) :: mix
151 : ! *************************************************************************
152 :
153 : ! Nullify internal pointers.
154 16160 : nullify(mix%i_rhor)
155 16160 : nullify(mix%i_vtrial)
156 16160 : nullify(mix%i_vresid)
157 16160 : nullify(mix%i_vrespc)
158 16160 : nullify(mix%f_fftgr)
159 16160 : nullify(mix%f_fftgr_sp)
160 16160 : nullify(mix%f_fftgr_trial_sp)
161 16160 : nullify(mix%f_fftgr_delta_i2)
162 16160 : nullify(mix%f_fftgr_delta_scale)
163 16160 : nullify(mix%f_atm)
164 16160 : nullify(mix%f_paw)
165 16160 : nullify(mix%f_extfpmd)
166 16160 : 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 8080 : 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 8080 : call init_(mix)
206 :
207 : ! Argument checkings.
208 8080 : 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 8080 : 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 8080 : & 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 8080 : errid = AB7_NO_ERROR
234 :
235 : ! Mandatory arguments.
236 8080 : mix%iscf = iscf
237 8080 : mix%kind = kind
238 8080 : mix%space = space
239 8080 : mix%nfft = nfft
240 8080 : mix%nspden = nspden
241 8080 : mix%n_pawmix = npawmix
242 :
243 : ! Optional arguments.
244 8080 : if (present(useprec)) mix%useprec = useprec
245 8080 : if (present(pulayhist_storage)) mix%pulayhist_storage = pulayhist_storage
246 8080 : 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 7925 : if (present(npulayit)) mix%n_pulayit = npulayit
300 7925 : mix%n_fftgr=2+2*mix%n_pulayit ; mix%n_index=1+mix%n_pulayit
301 7925 : if (.not. mix%useprec) mix%n_fftgr = 1+2*mix%n_pulayit
302 7925 : 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 24240 : ABI_MALLOC(mix%i_rhor,(mix%n_index))
315 71751 : 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 16160 : ABI_MALLOC(mix%i_vtrial,(mix%n_index))
319 71751 : 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 16160 : ABI_MALLOC(mix%i_vresid,(mix%n_index))
323 71751 : 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 16160 : ABI_MALLOC(mix%i_vrespc,(mix%n_index))
327 71751 : 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 63436 : do ii=1,mix%n_pulayit
360 63436 : mix%i_vtrial(ii)=2*ii-1 ; mix%i_vrespc(ii)=2*ii
361 : end do
362 7925 : mix%i_vrespc(mix%n_pulayit+1)=2*mix%n_pulayit+1
363 7925 : mix%i_vresid(1)=2*mix%n_pulayit+2
364 7925 : if (.not. mix%useprec) mix%i_vresid(1)=2
365 : end if ! iscf cases
366 :
367 8080 : 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 55346 : subroutine abi_mixing_copy_current_step(mix, arr_resid, errid, errmess, &
455 55456 : & 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 55346 : 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 55346 : 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 55346 : 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 55346 : 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 55346 : errid = AB7_NO_ERROR
502 :
503 55346 : if (mix%n_fftgr>0) then
504 55346 : 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 1272276394 : if (mix%i_vresid(1)>0) mix%f_fftgr(:,:,mix%i_vresid(1)) = arr_resid(:,:)
510 1240710524 : 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 55346 : if (mix%n_pawmix>0) then
514 14464 : if (present(arr_paw_resid).and.mix%i_vresid(1)>0) mix%f_paw(:, mix%i_vresid(1)) = arr_paw_resid(:)
515 14464 : 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 55346 : 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 55346 : if(associated(mix%f_extfpmd)) then
521 55346 : if(present(arr_extfpmd_resid).and.mix%i_vresid(1)>0) mix%f_extfpmd(mix%i_vresid(1))=arr_extfpmd_resid
522 55346 : if(present(arr_extfpmd_respc).and.mix%i_vrespc(1)>0) mix%f_extfpmd(mix%i_vrespc(1))=arr_extfpmd_respc
523 : endif
524 55346 : 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 110802 : 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 55346 : 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 55346 : istep_ = 1
566 55346 : if (present(istep)) istep_ = istep
567 :
568 : ! Allocate work array.
569 55346 : 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 7820 : 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 39095 : ABI_MALLOC(mix%f_fftgr,(mix%space * mix%nfft,mix%nspden,mix%n_fftgr))
577 : end if
578 2748274079 : mix%f_fftgr(:,:,:)=zero
579 7820 : 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 55346 : 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 55346 : 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 31280 : ABI_MALLOC(mix%f_paw,(mix%n_pawmix,mix%n_fftgr))
614 7820 : if (mix%n_pawmix > 0 .and. mix%n_fftgr>0) then
615 5376646 : mix%f_paw(:,:)=zero
616 1490 : 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 55346 : 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 31280 : ABI_MALLOC(mix%f_atm,(3,mix%n_atom,mix%n_fftgr))
628 : end if
629 :
630 55346 : if(.not.associated(mix%f_extfpmd)) then
631 7820 : if(mix%useextfpmd==1) then
632 21 : ABI_MALLOC(mix%f_extfpmd,(mix%n_fftgr))
633 119 : mix%f_extfpmd=zero
634 : else
635 7813 : ABI_MALLOC(mix%f_extfpmd,(0))
636 : endif
637 : endif
638 :
639 55346 : if(.not.associated(mix%f_rcpaw)) then
640 7820 : 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 7817 : ABI_MALLOC(mix%f_rcpaw,(0,0))
645 : endif
646 : endif
647 :
648 55346 : 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 55346 : 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 55346 : 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 55346 : 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 110692 : subroutine abi_mixing_eval(mix, arr, istep, nfftot, ucvol, &
738 : & mpi_comm, mpi_summarize, errid, errmess, &
739 54743 : & reset, isecur, pawarr, pawopt, response, etotal, potden, &
740 12922 : & 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 55346 : 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 14464 : 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 55346 : 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 55346 : errid = AB7_NO_ERROR
794 :
795 : ! Reset if requested
796 55346 : initialized = 1
797 55346 : if (present(reset)) then
798 55346 : if (reset) initialized = 0
799 : end if
800 :
801 : ! Miscellaneous
802 55346 : moveAtm = 0
803 55346 : if (mix%n_atom > 0) moveAtm = 1
804 55346 : isecur_ = 0
805 55346 : if (present(isecur)) isecur_ = isecur
806 55346 : comm_atom_ = xmpi_comm_self
807 55346 : if (present(comm_atom)) comm_atom_ = comm_atom
808 55346 : usepaw = 0
809 55346 : if (mix%n_pawmix > 0) usepaw = 1
810 55346 : pawoptmix_ = 0
811 55346 : if (present(pawopt)) pawoptmix_ = pawopt
812 55346 : response_ = 0
813 55346 : if (present(response)) response_ = response
814 55346 : pawarr_ => dum ; if (present(pawarr)) pawarr_ => pawarr
815 55346 : nelect_extfpmd_=zero;if(present(nelect_extfpmd)) nelect_extfpmd_=nelect_extfpmd
816 55346 : rcpawarr_ => dum0; if(present(rcpaw_arr)) rcpawarr_ => rcpaw_arr
817 :
818 : ! Do the mixing.
819 55346 : resnrm_ = 0.d0
820 55346 : if (mix%iscf == ABI_MIXING_NONE) then
821 31566473 : arr(:,:)=arr(:,:)+mix%f_fftgr(:,:,1)
822 54743 : 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 54696 : & mix%iscf == ABI_MIXING_ANDERSON_2 .or. &
830 : & mix%iscf == ABI_MIXING_PULAY) then
831 54510 : 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 54501 : 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 54501 : & 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 54510 : 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 55346 : if(present(nelect_extfpmd)) nelect_extfpmd=nelect_extfpmd_
894 55346 : if (present(resnrm)) resnrm = resnrm_
895 : nullify(rcpawarr_)
896 :
897 123011 : 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 8080 : 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 8080 : ABI_SFREE_PTR(mix%i_rhor)
929 8080 : ABI_SFREE_PTR(mix%i_vtrial)
930 8080 : ABI_SFREE_PTR(mix%i_vresid)
931 8080 : ABI_SFREE_PTR(mix%i_vrespc)
932 8080 : ABI_SFREE_PTR(mix%f_fftgr)
933 8080 : ABI_SFREE_PTR(mix%f_fftgr_sp)
934 8080 : ABI_SFREE_PTR(mix%f_fftgr_trial_sp)
935 8080 : ABI_SFREE_PTR(mix%f_fftgr_delta_i2)
936 8080 : ABI_SFREE_PTR(mix%f_fftgr_delta_scale)
937 8080 : ABI_SFREE_PTR(mix%f_paw)
938 8080 : ABI_SFREE_PTR(mix%f_atm)
939 8080 : ABI_SFREE_PTR(mix%f_extfpmd)
940 8080 : ABI_SFREE_PTR(mix%f_rcpaw)
941 :
942 8080 : call nullify_(mix)
943 :
944 8080 : 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 54501 : subroutine scfopt(cplex,f_fftgr,f_paw,iscf,istep,i_vrespc,i_vtrial,&
2012 : & mpicomm,mpi_summarize,nfft,npawmix,nspden,n_fftgr,&
2013 54501 : & n_index,opt_denpot,pawoptmix,usepaw,vpaw,vresid,vtrial,errid,errmess, &
2014 54501 : & useextfpmd,f_extfpmd,nelect_extfpmd,&
2015 54501 : & 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 54501 : 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 54501 : 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 54501 : errid = AB7_NO_ERROR
2060 :
2061 54501 : comm_atom_=xmpi_comm_self; if(present(comm_atom)) comm_atom_=comm_atom
2062 :
2063 54501 : i_vstore=i_vtrial(1)
2064 54501 : if (iscf==4) i_vstore=i_vtrial(2)
2065 54501 : if (iscf==7) then
2066 53990 : if (modulo(n_fftgr, 2) == 0 ) then
2067 53990 : npulay=(n_fftgr-2)/2
2068 : else
2069 0 : npulay=(n_fftgr-1)/2
2070 : end if
2071 53990 : 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 54501 : call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr)
2078 54501 : 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 54501 : vresid = resid_new(1)
2085 :
2086 : !_______________________________________________________________
2087 : !Here use only the preconditioning, or initialize the other algorithms
2088 :
2089 54501 : if (istep==1 .or. iscf==2) then
2090 8102 : write(message,'(2a)') ch10,' Simple mixing update:'
2091 8102 : call wrtout(std_out,message,'COLL')
2092 :
2093 8102 : write(message,*)' residual square of the potential: ',resid_new(1)
2094 8102 : call wrtout(std_out,message,'COLL')
2095 :
2096 : ! Store information for later use
2097 8102 : if (iscf==3.or.iscf==4) resid_old=resid_new(1)
2098 8102 : if (iscf==7) then
2099 7816 : amat(:,:)=zero
2100 7816 : amat(1,1)=resid_new(1)
2101 : end if
2102 :
2103 : ! Compute new vtrial (and new rhoij if PAW)
2104 179675876 : if (iscf/=2) f_fftgr(:,:,i_vstore)=vtrial(:,:)
2105 184891136 : vtrial(:,:)=vtrial(:,:)+f_fftgr(:,:,i_vrespc(1))
2106 8102 : if (usepaw==1) then
2107 381273 : if (iscf/=2) f_paw(:,i_vstore)=vpaw(:)
2108 382569 : vpaw(:)=vpaw(:)+f_paw(:,i_vrespc(1))
2109 : end if
2110 :
2111 8102 : 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 8102 : 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 46399 : 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 46197 : 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 46174 : else if(iscf==7)then
2284 :
2285 46174 : niter=min(istep,npulay+1)
2286 :
2287 46174 : write(message,'(2a,i2,a)') ch10,' Pulay update with ',niter-1,' previous iterations:'
2288 46174 : call wrtout(std_out,message,'COLL')
2289 :
2290 46174 : 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 46174 : if (istep>npulay+1) then
2300 83743 : do jj=1,niter-1
2301 597844 : do ii=1,niter-1
2302 587362 : amat(ii,jj)=amat(ii+1,jj+1)
2303 : end do
2304 : end do
2305 : end if
2306 46174 : 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 283782 : 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 237608 : & nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2319 237608 : if (usepaw==1.and.pawoptmix==1) amat(ii,niter)=amat(ii,niter)+amat_paw(ii)
2320 283782 : if (ii<niter) amat(niter,ii)=amat(ii,niter)
2321 : end do
2322 46174 : if (usepaw==1.and.pawoptmix==1)then
2323 143 : ABI_FREE(amat_paw)
2324 : end if
2325 :
2326 : ! Invert "A" matrix
2327 184696 : ABI_MALLOC(amatinv,(niter,niter))
2328 1741346 : amatinv(1:niter,1:niter)=amat(1:niter,1:niter)
2329 138522 : ABI_MALLOC(ipiv,(niter))
2330 138522 : ABI_MALLOC(rwork,(niter))
2331 46174 : call dgetrf(niter,niter,amatinv,niter,ipiv,ierr)
2332 46174 : call dgetri(niter,amatinv,niter,ipiv,rwork,niter,ierr)
2333 46174 : ABI_FREE(ipiv)
2334 46174 : ABI_FREE(rwork)
2335 :
2336 : ! Compute "alpha" factors
2337 138522 : ABI_MALLOC(alpha,(niter))
2338 46174 : det=zero
2339 283782 : alpha(:)=zero
2340 283782 : do ii=1,niter
2341 1741346 : do jj=1,niter
2342 1457564 : alpha(ii)=alpha(ii)+amatinv(jj,ii)
2343 1695172 : det=det+amatinv(jj,ii)
2344 : end do
2345 : end do
2346 283782 : alpha(:)=alpha(:)/det
2347 46174 : ABI_FREE(amatinv)
2348 46174 : 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 46174 : call wrtout(std_out,message,'COLL')
2350 :
2351 : ! Save latest trial potential and compute new trial potential
2352 104378 : do isp=1,nspden
2353 1050779846 : do ifft=1,cplex*nfft
2354 1050675468 : current=vtrial(ifft,isp)
2355 1050675468 : vtrial(ifft,isp)=alpha(niter)*(current+f_fftgr(ifft,isp,i_vrespc(1)))
2356 5783974768 : do ii=niter-1,1,-1
2357 : vtrial(ifft,isp)=vtrial(ifft,isp)+alpha(ii) &
2358 5783974768 : & *(f_fftgr(ifft,isp,i_vtrial(niter-ii))+f_fftgr(ifft,isp,i_vrespc(1+niter-ii)))
2359 : end do
2360 1050733672 : f_fftgr(ifft,isp,i_vstore)=current
2361 : end do
2362 : end do
2363 :
2364 : ! PAW: save latest rhoij and compute new rhoij
2365 3982395 : do index=1,npawmix
2366 3936221 : current=vpaw(index)
2367 3936221 : vpaw(index)=alpha(niter)*(current+f_paw(index,i_vrespc(1)))
2368 25175217 : do ii=niter-1,1,-1
2369 : vpaw(index)=vpaw(index)+alpha(ii) &
2370 25175217 : & *(f_paw(index,i_vtrial(niter-ii))+f_paw(index,i_vrespc(1+niter-ii)))
2371 : end do
2372 3982395 : f_paw(index,i_vstore)=current
2373 : end do
2374 :
2375 46174 : 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 46204 : 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 46204 : f_rcpaw(index,i_vstore)=current
2393 : end do
2394 :
2395 :
2396 46174 : ABI_FREE(alpha)
2397 : ! _______________________________________________________________
2398 : ! End of choice of optimization method
2399 : end if
2400 :
2401 : !Permute potential indices
2402 54501 : if (iscf==3) then
2403 216 : tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
2404 54285 : 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 54258 : else if (iscf==7) then
2408 53990 : tmp=i_vtrial( npulay)
2409 378132 : do ii= npulay,2,-1
2410 378132 : i_vtrial(ii)=i_vtrial(ii-1)
2411 : end do
2412 53990 : i_vtrial(1)=tmp
2413 53990 : tmp=i_vrespc(1+npulay)
2414 432122 : do ii=1+npulay,2,-1
2415 432122 : i_vrespc(ii)=i_vrespc(ii-1)
2416 : end do
2417 53990 : i_vrespc(1)=tmp
2418 : end if
2419 :
2420 54501 : 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 237856 : subroutine dotprodm_v(cplex,cpldot,dot,index1,index2,mpicomm,mpi_summarize,&
3043 237856 : & 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 237856 : if(cplex==1 .or. cpldot==1)then
3074 :
3075 475712 : do i1=1,mult1
3076 713568 : do i2=1,mult2
3077 237856 : ar=zero
3078 521602 : 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 5466280178 : do ifft=1,cplex*nfft
3081 5466042322 : ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1)
3082 : end do
3083 : end do
3084 237856 : dot(1,i1,i2)=ar
3085 475712 : if (nspden==4) then
3086 : ar=zero
3087 48390 : do ispden=3,4
3088 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
3089 320762374 : do ifft=1,cplex*nfft
3090 320746244 : ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1)
3091 : end do
3092 : end do
3093 16130 : if (opt_storage==0) then
3094 9392 : if (cplex==1) then
3095 8844 : 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 237856 : if (mpi_summarize) then
3146 31776 : call timab(48,1,tsec)
3147 31776 : call xmpi_sum(dot,mpicomm ,ierr)
3148 31776 : call timab(48,2,tsec)
3149 : end if
3150 :
3151 237856 : if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero
3152 :
3153 237856 : 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 54687 : 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 109374 : do ii=1,mult
3523 54687 : ar=zero
3524 117785 : do ispden=1,min(nspden,2)
3525 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar)
3526 1188997400 : do ifft=1,cplex*nfft
3527 1188942713 : ar=ar + potarr(ifft,ispden,index+ii-1)**2
3528 : end do
3529 : end do
3530 54687 : norm2(ii)=ar
3531 109374 : if (nspden==4) then
3532 : ar=zero
3533 7524 : do ispden=3,4
3534 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar)
3535 51154280 : do ifft=1,cplex*nfft
3536 51151772 : ar=ar + potarr(ifft,ispden,index+ii-1)**2
3537 : end do
3538 : end do
3539 2508 : if (opt_storage==0) then
3540 1517 : if (cplex==1) then
3541 1412 : 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 54687 : if (mpi_summarize) then
3553 6780 : call timab(48,1,tsec)
3554 6780 : call xmpi_sum(norm2,mpicomm ,ierr)
3555 6780 : call timab(48,2,tsec)
3556 : end if
3557 :
3558 54687 : 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 : !!***
|
