TABLE OF CONTENTS


ABINIT/scfopt [ Functions ]

[ Top ] [ Functions ]

NAME

 scfopt

FUNCTION

 Compute the next vtrial of the SCF cycle.
 Possible algorithms are : simple mixing, Anderson (order 1 or 2), Pulay

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (XG,GMR)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
  iscf= 2 => simple mixing
      = 3,4 => Anderson mixing
      = 7 => Pulay mixing
  istep= number of the step in the SCF cycle
  mpicomm=the mpi communicator used for the summation
  comm_atom=the mpi communicator over atoms ; PAW only (optional argument)
  mpi_summarize=set it to .true. if parallelisation is done over FFT
  nfft=(effective) number of FFT grid points (for this processor)
  npawmix=-PAW only- number of spherical part elements to be mixed
  nspden=number of spin-density components
  n_fftgr=third dimension of the array f_fftgr
  n_index=dimension for indices of potential/density (see ivrespc, i_vtrial...)
  opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential
              1 vtrial (and also f_fftgr) actually contains the trial density
  pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part
  usepaw= 0 for non paw calculation; =1 for paw calculation

OUTPUT

  (see side effects)

SIDE EFFECTS

  vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave
     the input preconditioned residual potential
     at output, it is the new trial potential .
  f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid :
   The input vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(1)).
   The old vtrial is transferred, at output,in f_fftgr(:,:,i_vtrial(2)).
   The input preconditioned residual potential is in f_fftgr(:,:,i_vrespc(1))
   Two input old preconditioned residual potentials in f_fftgr(:,:,i_vrespc(2)) and f_fftgr(:,:,i_vrespc(3))
    Before output a permutation of i_vrespc(1), i_vrespc(2) and i_vrespc(3) occurs, without
    actually copying all the data (change of pointer).
  i_vrespc(n_index)=index of the preconditioned residual potentials (present and past) in the array f_fftgr
  i_vtrial(n_index)  =indices of the potential (present and past) in the array f_fftgr
  ==== if usepaw==1
    f_paw(npawmix,n_fftgr*mffmem*usepaw)=different functions used for PAW
                                           (same as f_fftgr but for spherical part)
    vpaw(npawmix*usepaw)=at input, the aug. occupancies (rhoij) that gave
                               the input preconditioned residual potential
                           at output, it is the new aug. occupancies.

PARENTS

      m_ab7_mixing

CHILDREN

      dgetrf,dgetri,dotprodm_v,sqnormm_v,wrtout,xmpi_sum

SOURCE

 68 #if defined HAVE_CONFIG_H
 69 #include "config.h"
 70 #endif
 71 
 72 #include "abi_common.h"
 73 
 74 subroutine scfopt(cplex,f_fftgr,f_paw,iscf,istep,i_vrespc,i_vtrial,&
 75 & mpicomm,mpi_summarize,nfft,npawmix,nspden,n_fftgr,&
 76 & n_index,opt_denpot,pawoptmix,usepaw,vpaw,vresid,vtrial,errid,errmess, &
 77 & comm_atom) ! optional
 78 
 79  use defs_basis
 80  use m_linalg_interfaces
 81  use m_profiling_abi
 82  use m_xmpi
 83 
 84 !This section has been created automatically by the script Abilint (TD).
 85 !Do not modify the following lines by hand.
 86 #undef ABI_FUNC
 87 #define ABI_FUNC 'scfopt'
 88  use interfaces_14_hidewrite
 89  use interfaces_56_mixing, except_this_one => scfopt
 90 !End of the abilint section
 91 
 92  implicit none
 93 
 94 !Arguments ------------------------------------
 95 !scalars
 96  integer,intent(in) :: cplex,iscf,istep,n_fftgr,n_index,nfft
 97  integer,intent(in) :: npawmix,nspden,opt_denpot,pawoptmix,usepaw,mpicomm
 98  integer, intent(in),optional :: comm_atom
 99  integer,intent(out) :: errid
100  character(len = 500), intent(out) :: errmess
101  logical, intent(in) :: mpi_summarize
102  real(dp), intent(out) :: vresid
103 !arrays
104  integer,intent(inout) :: i_vrespc(n_index),i_vtrial(n_index)
105  real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr)
106  real(dp),intent(inout) :: f_paw(npawmix,n_fftgr*usepaw),vpaw(npawmix*usepaw)
107  real(dp),intent(inout) :: vtrial(cplex*nfft,nspden)
108 !Local variables-------------------------------
109 !scalars
110  integer,parameter :: npulaymax=50
111  integer :: i_vstore,ierr,ifft,ii,index,isp,jj,comm_atom_,niter,npulay,tmp
112  real(dp),save :: prod_resid_old,resid_old,resid_old2
113  real(dp) :: aa1,aa2,bb,cc1,cc2,current,det,lambda,lambda2,resid_best
114  character(len=500) :: message
115 !arrays
116  integer,allocatable :: ipiv(:)
117  real(dp),save :: amat(npulaymax+1,npulaymax+1)
118  real(dp) :: mpibuff(2),prod_resid(1),prod_resid2(1),resid_new(1)
119  real(dp),allocatable :: alpha(:),amatinv(:,:),amat_paw(:),rwork(:)
120 
121 ! *************************************************************************
122 
123 !DEBUG
124 !write(std_out,*)' scfopt : enter ; istep,iscf ',istep,iscf
125 !ENDDEBUG
126 
127  errid = AB7_NO_ERROR
128 
129  comm_atom_=xmpi_comm_self; if(present(comm_atom)) comm_atom_=comm_atom
130 
131  i_vstore=i_vtrial(1)
132  if (iscf==4) i_vstore=i_vtrial(2)
133  if (iscf==7) then
134    if (modulo(n_fftgr, 2) == 0 ) then
135      npulay=(n_fftgr-2)/2
136    else
137      npulay=(n_fftgr-1)/2
138    end if
139    i_vstore=i_vtrial(npulay)
140  else
141    npulay=0
142  end if
143 
144 !Compute the new residual resid_new, from f_fftgr/f_paw(:,:,i_vrespc(1))
145  call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr)
146  if (usepaw==1.and.pawoptmix==1) then
147    do index=1,npawmix
148      resid_new(1)=resid_new(1)+f_paw(index,i_vrespc(1))**2
149    end do
150    call xmpi_sum(resid_new(1),comm_atom_,ierr)
151  end if
152  vresid = resid_new(1)
153 
154 !_______________________________________________________________
155 !Here use only the preconditioning, or initialize the other algorithms
156 
157  if(istep==1 .or. iscf==2)then
158 
159    write(message,'(2a)') ch10,'Simple mixing update:'
160    call wrtout(std_out,message,'COLL')
161 
162    write(message,*)' residual square of the potential :',resid_new(1)
163    call wrtout(std_out,message,'COLL')
164 
165 !  Store information for later use
166    if (iscf==3.or.iscf==4) resid_old=resid_new(1)
167    if (iscf==7) then
168      amat(:,:)=zero
169      amat(1,1)=resid_new(1)
170    end if
171 
172 !  Compute new vtrial (and new rhoij if PAW)
173    if (iscf/=2) f_fftgr(:,:,i_vstore)=vtrial(:,:)
174    vtrial(:,:)=vtrial(:,:)+f_fftgr(:,:,i_vrespc(1))
175    if (usepaw==1) then
176      if (iscf/=2) f_paw(:,i_vstore)=vpaw(:)
177      vpaw(:)=vpaw(:)+f_paw(:,i_vrespc(1))
178    end if
179 
180 !  _______________________________________________________________
181 !  Here Anderson algorithm using one previous iteration
182  else if((istep==2 .or. iscf==3).and.iscf/=7)then
183 
184    write(message,'(2a)') ch10,'Anderson update:'
185    call wrtout(std_out,message,'COLL')
186 
187    write(message,*)' residual square of the potential: ',resid_new(1)
188    call wrtout(std_out,message,'COLL')
189 
190 !  Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2))
191    call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,&
192 &   nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
193    if (usepaw==1.and.pawoptmix==1) then
194      do index=1,npawmix
195        prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2))
196      end do
197      call xmpi_sum(prod_resid(1),comm_atom_,ierr)
198    end if
199 
200 !  Compute mixing factor
201    lambda=(resid_new(1)-prod_resid(1))/(resid_new(1)+resid_old-2*prod_resid(1))
202    write(message,*)' mixing of old trial potential    :',lambda
203    call wrtout(std_out,message,'COLL')
204 
205 !  Evaluate best residual square on the line
206    resid_best=(1.0_dp-lambda)*(1.0_dp-lambda)*resid_new(1)&
207 &   +(1.0_dp-lambda)*lambda        *2*prod_resid(1)&
208 &   +lambda        *lambda        *resid_old
209    write(message,*)' predicted best residual square on the line: ',resid_best
210    call wrtout(std_out,message,'COLL')
211 
212 !  Store information for later use
213    if (iscf==4) then
214      prod_resid_old=prod_resid(1)
215      resid_old2=resid_old
216    end if
217    resid_old=resid_new(1)
218 
219 !  Save latest trial potential and compute new trial potential
220    do isp=1,nspden
221      do ifft=1,cplex*nfft
222        current=vtrial(ifft,isp)
223        vtrial(ifft,isp)=(one-lambda)*(current                      +f_fftgr(ifft,isp,i_vrespc(1)))&
224 &       +lambda      *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))
225        f_fftgr(ifft,isp,i_vstore)=current
226      end do
227    end do
228 
229 !  PAW: save latest rhoij and compute new rhoij
230    do index=1,npawmix
231      current=vpaw(index)
232      vpaw(index)=(one-lambda)*(current                 +f_paw(index,i_vrespc(1)))&
233 &     +lambda      *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))
234      f_paw(index,i_vstore)=current
235    end do
236 
237 !  _______________________________________________________________
238 !  Here Anderson algorithm using two previous iterations
239  else if(iscf==4.and.iscf/=7)then
240 
241    write(message,'(2a)') ch10,'Anderson (order 2) update:'
242    call wrtout(std_out,message,'COLL')
243 
244    write(message,*)' residual square of the potential :',resid_new(1)
245    call wrtout(std_out,message,'COLL')
246 
247 !  Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2))
248    call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,&
249 &   nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
250    if (usepaw==1.and.pawoptmix==1) then
251      do index=1,npawmix
252        prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2))
253      end do
254    end if
255 
256 !  Compute prod_resid2 from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(3))
257    call dotprodm_v(cplex,1,prod_resid2,i_vrespc(1),i_vrespc(3),mpicomm,mpi_summarize,1,1,&
258 &   nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
259    if (usepaw==1.and.pawoptmix==1) then
260      do index=1,npawmix
261        prod_resid2(1)=prod_resid2(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(3))
262      end do
263 !    MPI reduction
264      mpibuff(1)=prod_resid(1);mpibuff(2)=prod_resid2(1)
265      call xmpi_sum(mpibuff,comm_atom_,ierr)
266      prod_resid(1)=mpibuff(1);prod_resid2(1)=mpibuff(2)
267    end if
268 
269 !  Compute mixing factors
270    aa1=resid_new(1)+resid_old -two*prod_resid (1)
271    aa2=resid_new(1)+resid_old2-two*prod_resid2(1)
272    bb =resid_new(1)+prod_resid_old-prod_resid(1)-prod_resid2(1)
273    cc1=resid_new(1)-prod_resid (1)
274    cc2=resid_new(1)-prod_resid2(1)
275    det=aa1*aa2-bb*bb
276    lambda =(aa2*cc1-bb*cc2)/det
277    lambda2=(aa1*cc2-bb*cc1)/det
278    write(message,*)' mixing of old trial potentials   :',lambda,lambda2
279    call wrtout(std_out,message,'COLL')
280 
281 !  Store information for later use
282    prod_resid_old=prod_resid(1)
283    resid_old2=resid_old
284    resid_old=resid_new(1)
285 
286 !  Save latest trial potential and compute new trial potential
287    do isp=1,nspden
288      do ifft=1,cplex*nfft
289        current=vtrial(ifft,isp)
290        vtrial(ifft,isp)=&
291 &       (one-lambda-lambda2)*(current                      +f_fftgr(ifft,isp,i_vrespc(1)))&
292 &       +lambda             *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))&
293 &       +lambda2            *(f_fftgr(ifft,isp,i_vtrial(2))+f_fftgr(ifft,isp,i_vrespc(3)))
294        f_fftgr(ifft,isp,i_vstore)=current
295      end do
296    end do
297 
298 !  PAW: save latest rhoij and compute new rhoij
299    do index=1,npawmix
300      current=vpaw(index)
301      vpaw(index)=&
302 &     (one-lambda-lambda2)*(current                 +f_paw(index,i_vrespc(1)))&
303 &     +lambda             *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))&
304 &     +lambda2            *(f_paw(index,i_vtrial(2))+f_paw(index,i_vrespc(3)))
305      f_paw(index,i_vstore)=current
306    end do
307 
308 !  _______________________________________________________________
309 !  Here Pulay algorithm
310  else if(iscf==7)then
311 
312    niter=min(istep,npulay+1)
313 
314    write(message,'(2a,i2,a)') ch10,' Pulay update with ',niter-1,' previous iterations:'
315    call wrtout(std_out,message,'COLL')
316 
317    if (npulay>npulaymax) then
318      errid = AB7_ERROR_MIXING_CONVERGENCE
319      write(errmess, '(4a)' ) ch10,&
320 &     ' scfopt : ERROR - ',ch10,&
321 &     '  Too much iterations required for Pulay algorithm (<50) !'
322      return
323    end if
324 
325 !  Compute "A" matrix
326    if (istep>npulay+1) then
327      do jj=1,niter-1
328        do ii=1,niter-1
329          amat(ii,jj)=amat(ii+1,jj+1)
330        end do
331      end do
332    end if
333    if (usepaw==1.and.pawoptmix==1) then
334      ABI_ALLOCATE(amat_paw,(niter))
335      amat_paw(:)=zero
336      do ii=1,niter
337        do index=1,npawmix
338          amat_paw(ii)=amat_paw(ii)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(1+niter-ii))
339        end do
340      end do
341      call xmpi_sum(amat_paw,comm_atom_,ierr)
342    end if
343    do ii=1,niter
344      call dotprodm_v(cplex,1,amat(ii,niter),i_vrespc(1),i_vrespc(1+niter-ii),mpicomm,mpi_summarize,1,1,&
345 &     nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
346      if (usepaw==1.and.pawoptmix==1) amat(ii,niter)=amat(ii,niter)+amat_paw(ii)
347      if (ii<niter) amat(niter,ii)=amat(ii,niter)
348    end do
349    if (usepaw==1.and.pawoptmix==1)then
350      ABI_DEALLOCATE(amat_paw)
351    end if
352 
353 !  Invert "A" matrix
354    ABI_ALLOCATE(amatinv,(niter,niter))
355    amatinv(1:niter,1:niter)=amat(1:niter,1:niter)
356    ABI_ALLOCATE(ipiv,(niter))
357    ABI_ALLOCATE(rwork,(niter))
358    call dgetrf(niter,niter,amatinv,niter,ipiv,ierr)
359    call dgetri(niter,amatinv,niter,ipiv,rwork,niter,ierr)
360    ABI_DEALLOCATE(ipiv)
361    ABI_DEALLOCATE(rwork)
362 
363 !  Compute "alpha" factors
364    ABI_ALLOCATE(alpha,(niter))
365    alpha=zero
366    det=zero
367    do ii=1,niter
368      do jj=1,niter
369        alpha(ii)=alpha(ii)+amatinv(jj,ii)
370        det=det+amatinv(jj,ii)
371      end do
372    end do
373    alpha(:)=alpha(:)/det
374    ABI_DEALLOCATE(amatinv)
375    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)
376    call wrtout(std_out,message,'COLL')
377 
378 !  Save latest trial potential and compute new trial potential
379    do isp=1,nspden
380      do ifft=1,cplex*nfft
381        current=vtrial(ifft,isp)
382        vtrial(ifft,isp)=alpha(niter)*(current+f_fftgr(ifft,isp,i_vrespc(1)))
383        do ii=niter-1,1,-1
384          vtrial(ifft,isp)=vtrial(ifft,isp)+alpha(ii) &
385 &         *(f_fftgr(ifft,isp,i_vtrial(niter-ii))+f_fftgr(ifft,isp,i_vrespc(1+niter-ii)))
386        end do
387        f_fftgr(ifft,isp,i_vstore)=current
388      end do
389    end do
390 
391 !  PAW: save latest rhoij and compute new rhoij
392    do index=1,npawmix
393      current=vpaw(index)
394      vpaw(index)=alpha(niter)*(current+f_paw(index,i_vrespc(1)))
395      do ii=niter-1,1,-1
396        vpaw(index)=vpaw(index)+alpha(ii) &
397 &       *(f_paw(index,i_vtrial(niter-ii))+f_paw(index,i_vrespc(1+niter-ii)))
398      end do
399      f_paw(index,i_vstore)=current
400    end do
401 
402    ABI_DEALLOCATE(alpha)
403 !  _______________________________________________________________
404 !  End of choice of optimization method
405  end if
406 
407 !Permute potential indices
408  if (iscf==3) then
409    tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
410  else if (iscf==4) then
411    tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
412    tmp=i_vtrial(2) ; i_vtrial(2)=i_vtrial(1) ; i_vtrial(1)=tmp
413  else if (iscf==7) then
414    tmp=i_vtrial(  npulay)
415    do ii=  npulay,2,-1
416      i_vtrial(ii)=i_vtrial(ii-1)
417    end do
418    i_vtrial(1)=tmp
419    tmp=i_vrespc(1+npulay)
420    do ii=1+npulay,2,-1
421      i_vrespc(ii)=i_vrespc(ii-1)
422    end do
423    i_vrespc(1)=tmp
424  end if
425 
426 end subroutine scfopt