TABLE OF CONTENTS


ABINIT/aprxdr [ Functions ]

[ Top ] [ Functions ]

NAME

 aprxdr

FUNCTION

 Compute the approximative derivatives of the energy at different
 points along the line search, thanks to a finite-difference formula.
 This formula is the projection along the line search of the
 Eq.(11) in PRB54, 4383 (1996) [[cite:Gonze1996]].

INPUTS

 cplex: if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
 choice= if==3, compute dedv_new, dedv_old, and dedv_mix,
 if/=3, compute only dedv_new and dedv_old.
 i_vresid and i_rhor, see the next lines.
 f_fftgr(nfft,nspden,n_fftgr)=different functions defined on the fft grid :
 The last residual potential is in f_fftgr(:,:,i_vresid(1)).
 The old  residual potential is in f_fftgr(:,:,i_vresid(2)).
 The previous old residual potential is in f_fftgr(:,:,i_vresid(3)).
 (needed only when choice==3)
 The old  density is in f_fftgr(:,:,i_rhor2).
 f_atm(3,natom,n_fftgr)=different functions defined for each atom :
 The last HF force is in f_atm(:,:,i_vresid(1)).
 The old  HF force is in f_fftgr(:,:,i_vresid(2)).
 The previous old HF force is in f_fftgr(:,:,i_vresid(3)).
 (needed only when choice==3)
 The old atomic positions are in f_atm(:,:,i_rhor2)
 moved_atm_inside: if==1, the atoms are allowed to move.
 mpicomm=the mpi communicator used for the summation
 mpi_summarize=set it to .true. if parallelisation is done over FFT
 natom=number of atoms in unit cell
 nfft=(effective) number of FFT grid points (for this processor)
 nfftot=total number of FFT grid points
 nspden=number of spin-density components
 rhor(nfft,nspden)=actual density
 xred(3,natom)=reduced atomic coordinates

OUTPUT

 dedv_mix=approximate derivative from previous old residual
 dedv_new=approximate derivative from new residual
 dedv_old=approximate derivative from old residual (output only when choice==3)

NOTES

 Should be OpenMP parallelized

PARENTS

      scfcge

CHILDREN

      dotprodm_vn

SOURCE

3183 subroutine aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,&
3184 &  f_atm,f_fftgr,i_rhor2,i_vresid,moved_atm_inside,&
3185 &  mpicomm,mpi_summarize,natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred)
3186 
3187 
3188 !This section has been created automatically by the script Abilint (TD).
3189 !Do not modify the following lines by hand.
3190 #undef ABI_FUNC
3191 #define ABI_FUNC 'aprxdr'
3192 !End of the abilint section
3193 
3194  implicit none
3195 
3196 !Arguments ------------------------------------
3197 !scalars
3198  integer,intent(in) :: choice,cplex,i_rhor2,moved_atm_inside,n_fftgr,natom,nfft
3199  integer,intent(in) :: mpicomm,nfftot,nspden
3200  logical, intent(in) :: mpi_summarize
3201  real(dp),intent(in) :: ucvol
3202  real(dp),intent(out) :: dedv_mix,dedv_new,dedv_old
3203 !arrays
3204  integer,intent(in) :: i_vresid(3)
3205  real(dp),intent(in) :: f_atm(3,natom,n_fftgr)
3206  real(dp),intent(in) :: f_fftgr(cplex*nfft,nspden,n_fftgr)
3207  real(dp),intent(in) :: rhor(cplex*nfft,nspden),xred(3,natom)
3208 
3209 !Local variables-------------------------------
3210 !scalars
3211  integer :: iatom,idir
3212 !arrays
3213  real(dp) :: dedv_temp(1)
3214  real(dp),allocatable :: ddens(:,:,:)
3215 
3216 ! *************************************************************************
3217 
3218  ABI_ALLOCATE(ddens,(cplex*nfft,nspden,1))
3219 
3220 !Compute approximative derivative of the energy
3221 !with respect to change of potential
3222 
3223  ddens(:,:,1)=rhor(:,:)-f_fftgr(:,:,i_rhor2)
3224 
3225 !call dotprod_vn(cplex,1,ddens,dedv_old,nfft,nfftot,nspden,1,vresid,ucvol)
3226 !Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(2))
3227  call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(2),mpicomm,mpi_summarize,1,1,1,&
3228 & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol)
3229  dedv_old = dedv_temp(1)
3230 
3231 !Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(1))
3232  call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(1),mpicomm,mpi_summarize,1,1,1,&
3233 & nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol)
3234  dedv_new= dedv_temp(1)
3235 
3236  if(choice==3)then
3237 !  Dot product ddens(:,:,1) f_fftgr(:,:,i_vresid(3))
3238    call dotprodm_vn(cplex,1,ddens,dedv_temp,1,i_vresid(3),mpicomm,mpi_summarize,1,1,1,&
3239 &   nfft,nfftot,n_fftgr,nspden,f_fftgr,ucvol)
3240    dedv_mix = dedv_temp(1)
3241  end if
3242 
3243  ABI_DEALLOCATE(ddens)
3244 
3245 !-------------------------------------------------------
3246 
3247 !Now, take care of eventual atomic displacements
3248 
3249  if(moved_atm_inside==1)then
3250    do idir=1,3
3251      do iatom=1,natom
3252        dedv_new=dedv_new+&
3253 &       f_atm(idir,iatom,i_vresid(1))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2))
3254        dedv_old=dedv_old+&
3255 &       f_atm(idir,iatom,i_vresid(2))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2))
3256        if(choice==3) dedv_mix=dedv_mix+&
3257 &       f_atm(idir,iatom,i_vresid(3))*(xred(idir,iatom)-f_atm(idir,iatom,i_rhor2))
3258      end do
3259    end do
3260  end if
3261 
3262 end subroutine aprxdr

ABINIT/dotprodm_v [ Functions ]

[ Top ] [ Functions ]

NAME

 dotprodm_v

FUNCTION

 For two sets of potentials,
 compute dot product of each pair of two potentials (integral over FFT grid), to obtain
 a series of square residual-like quantity (so the sum of product of values
 is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume).
 Take into account the spin components of the potentials (nspden),
 and sum over them.
 Need the index of the first pair of potentials to be treated, in each array
 of potentials, and the number of potentials to be treated.
 Might be used to compute just one square of norm, in
 a big array, such as to avoid copying a potential from a big array
 to a temporary place.

INPUTS

  cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
  cpldot=if 1, the dot array is real, if 2, the dot array is complex
  index1=index of the first potential to be treated in the potarr1 array
  index2=index of the first potential to be treated in the potarr2 array
  mpicomm=the mpi communicator used for the summation
  mpi_summarize=set it to .true. if parallelisation is done over FFT
  mult1=number of potentials to be treated in the first set
  mult2=number of potentials to be treated in the second set
  nfft= (effective) number of FFT grid points (for this processor)
  npot1= third dimension of the potarr1 array
  npot2= third dimension of the potarr2 array
  nspden=number of spin-density components
  opt_storage: 0, if potentials are stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn]
               1, if potentials are stored as V, B_x, B_y, Bz  (B=magn. field)
  potarr1(cplex*nfft,nspden,npot)=first array of real space potentials on FFT grid
    (if cplex=2 and cpldot=2, potarr1 is the array that will be complex conjugated)
  potarr2(cplex*nfft,nspden,npot)=second array of real space potentials on FFT grid

OUTPUT

  dot(cpldot,mult1,mult2)= series of values of the dot product

SIDE EFFECTS

NOTES

  Concerning storage when nspden=4:
   cplex=1:
     opt_storage=0: V are stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
     opt_storage=1: V are stored as : V, B_x, B_y, B_z               (real)
   cplex=2:
     opt_storage=0: V are stored as : V^11, V^22, V^12, i.V^21 (complex)
     opt_storage=1: V are stored as : V, B_x, B_y, B_z         (complex)

PARENTS

      scfopt

CHILDREN

      timab,xmpi_sum

SOURCE

2572 subroutine dotprodm_v(cplex,cpldot,dot,index1,index2,mpicomm,mpi_summarize,&
2573 &   mult1,mult2,nfft,npot1,npot2,nspden,opt_storage,potarr1,potarr2)
2574 
2575 
2576 !This section has been created automatically by the script Abilint (TD).
2577 !Do not modify the following lines by hand.
2578 #undef ABI_FUNC
2579 #define ABI_FUNC 'dotprodm_v'
2580 !End of the abilint section
2581 
2582  implicit none
2583 
2584 !Arguments ------------------------------------
2585 !scalars
2586  integer,intent(in) :: cpldot,cplex,index1,index2,mult1,mult2,nfft,npot1,npot2
2587  integer,intent(in) :: nspden,opt_storage,mpicomm
2588  logical, intent(in) :: mpi_summarize
2589 !arrays
2590  real(dp),intent(in) :: potarr1(cplex*nfft,nspden,npot1)
2591  real(dp),intent(in) :: potarr2(cplex*nfft,nspden,npot2)
2592  real(dp),intent(out) :: dot(cpldot,mult1,mult2)
2593 
2594 !Local variables-------------------------------
2595 !scalars
2596  integer :: i1,i2,ierr,ifft,ispden
2597  real(dp) :: ai,ar
2598 !arrays
2599  real(dp) :: tsec(2)
2600 
2601 ! *************************************************************************
2602 
2603 !Real or complex inputs are coded
2604  DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
2605 
2606 !Real or complex outputs are coded
2607  DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot")
2608  DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
2609  DBG_CHECK( npot1-index1-mult1 >= -1,"npot1-index1-mult1")
2610  DBG_CHECK( npot2-index2-mult2 >= -1,"npot2-index2-mult2")
2611 
2612  if(cplex==1 .or. cpldot==1)then
2613 
2614    do i1=1,mult1
2615      do i2=1,mult2
2616        ar=zero
2617        do ispden=1,min(nspden,2)
2618 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
2619          do ifft=1,cplex*nfft
2620            ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1)
2621          end do
2622        end do
2623        dot(1,i1,i2)=ar
2624        if (nspden==4) then
2625          ar=zero
2626          do ispden=3,4
2627 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar)
2628            do ifft=1,cplex*nfft
2629              ar=ar + potarr1(ifft,ispden,index1+i1-1)*potarr2(ifft,ispden,index2+i2-1)
2630            end do
2631          end do
2632          if (opt_storage==0) then
2633            if (cplex==1) then
2634              dot(1,i1,i2)=dot(1,i1,i2)+two*ar
2635            else
2636              dot(1,i1,i2)=dot(1,i1,i2)+ar
2637            end if
2638          else
2639            dot(1,i1,i2)=half*(dot(1,i1,i2)+ar)
2640          end if
2641        end if
2642      end do
2643    end do
2644 
2645  else ! if (cplex==2 .and. cpldot==2)
2646 
2647    do i1=1,mult1
2648      do i2=1,mult2
2649        ar=zero ; ai=zero
2650        do ispden=1,min(nspden,2)
2651 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
2652          do ifft=1,nfft
2653            ar=ar + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) &
2654 &           + potarr1(2*ifft  ,ispden,index1+i1-1)*potarr2(2*ifft  ,ispden,index2+i2-1)
2655            ai=ai + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft  ,ispden,index2+i2-1) &
2656 &           - potarr1(2*ifft  ,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1)
2657          end do
2658        end do
2659        dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai
2660        if (nspden==4) then
2661          ar=zero
2662          do ispden=3,4
2663 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,i1,i2,index1,index2,ispden,nfft,potarr1,potarr2) REDUCTION(+:ar,ai)
2664            do ifft=1,nfft
2665              ar=ar + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1) &
2666 &             + potarr1(2*ifft  ,ispden,index1+i1-1)*potarr2(2*ifft  ,ispden,index2+i2-1)
2667              ai=ai + potarr1(2*ifft-1,ispden,index1+i1-1)*potarr2(2*ifft  ,ispden,index2+i2-1) &
2668 &             - potarr1(2*ifft  ,ispden,index1+i1-1)*potarr2(2*ifft-1,ispden,index2+i2-1)
2669            end do
2670          end do
2671          if (opt_storage==0) then
2672            dot(1,i1,i2)=dot(1,i1,i2)+ar
2673            dot(2,i1,i2)=dot(2,i1,i2)+ai
2674          else
2675            dot(1,i1,i2)=half*(dot(1,i1,i2)+ar)
2676            dot(2,i1,i2)=half*(dot(2,i1,i2)+ai)
2677          end if
2678        end if
2679      end do
2680    end do
2681  end if
2682 
2683 !XG030513 : MPIWF reduction (addition) on dot is needed here
2684  if (mpi_summarize) then
2685    call timab(48,1,tsec)
2686    call xmpi_sum(dot,mpicomm ,ierr)
2687    call timab(48,2,tsec)
2688  end if
2689 
2690  if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero
2691 
2692 end subroutine dotprodm_v

ABINIT/dotprodm_vn [ Functions ]

[ Top ] [ Functions ]

NAME

 dotprodm_vn

FUNCTION

 For a set of densities and a set of potentials,
 compute the dot product (integral over FFT grid) of each pair, to obtain
 a series of energy-like quantity (so the usual dotproduct is divided
 by the number of FFT points, and multiplied by the primitive cell volume).
 Take into account the spin components of the density and potentials (nspden),
 and sum correctly over them. Note that the storage of densities and
 potentials is different : for potential, one stores the matrix components,
 while for the density, one stores the trace, and then, either the
 spin-polarisation (if nspden=2), or the magnetization vector (if nspden=4).
 Need the index of the first density/potential pair to be treated, in each array,
 and the number of pairs to be treated.
 Might be used to compute just one dot product, in
 a big array, such as to avoid copying the density and potential from a big array
 to a temporary place.

INPUTS

  cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
  cpldot=if 1, the dot array is real, if 2, the dot array is complex (not coded yet for nspden=4)
  denarr(cplex*nfft,nspden,nden)=real space density on FFT grid
  id=index of the first density to be treated in the denarr array
  ip=index of the first potential to be treated in the potarr array
  mpicomm=the mpi communicator used for the summation
  mpi_summarize=set it to .true. if parallelisation is done over FFT
  multd=number of densities to be treated
  multp=number of potentials to be treated
  nden=third dimension of the denarr array
  nfft= (effective) number of FFT grid points (for this processor)
  nfftot= total number of FFT grid points
  npot=third dimension of the potarr array
  nspden=number of spin-density components
  potarr(cplex*nfft,nspden,npot)=real space potential on FFT grid
                 (will be complex conjugated if cplex=2 and cpldot=2)
  ucvol=unit cell volume (Bohr**3)

OUTPUT

  dot(cpldot,multp,multd)= series of values of the dot product potential/density

SIDE EFFECTS

NOTES

  Concerning storage when nspden=4:
   cplex=1:
     V are stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
     N are stored as : n, m_x, m_y, m_z               (real)
   cplex=2:
     V are stored as : V^11, V^22, V^12, i.V^21 (complex)
     N are stored as : n, m_x, m_y, mZ          (complex)

PARENTS

      aprxdr

CHILDREN

      timab,xmpi_sum

SOURCE

2756 subroutine dotprodm_vn(cplex,cpldot,denarr,dot,id,ip,mpicomm, mpi_summarize,multd,multp,&
2757 & nden,nfft,nfftot,npot,nspden,potarr,ucvol)
2758 
2759 
2760 !This section has been created automatically by the script Abilint (TD).
2761 !Do not modify the following lines by hand.
2762 #undef ABI_FUNC
2763 #define ABI_FUNC 'dotprodm_vn'
2764 !End of the abilint section
2765 
2766  implicit none
2767 
2768 !Arguments ------------------------------------
2769 !scalars
2770  integer,intent(in) :: cpldot,cplex,id,ip,multd,multp,nden,nfft,nfftot,npot
2771  integer,intent(in) :: nspden,mpicomm
2772  logical, intent(in) :: mpi_summarize
2773  real(dp),intent(in) :: ucvol
2774 !arrays
2775  real(dp),intent(in) :: denarr(cplex*nfft,nspden,nden)
2776  real(dp),intent(in) :: potarr(cplex*nfft,nspden,npot)
2777  real(dp),intent(out) :: dot(cpldot,multp,multd)
2778 
2779 !Local variables-------------------------------
2780 !scalars
2781  integer :: i1,i2,ierr,ir,jr
2782  real(dp) :: ai,ar,dim11,dim12,dim21,dim22,dim_dn,dim_up,dre11,dre12,dre21
2783  real(dp) :: dre22,dre_dn,dre_up,factor,pim11,pim12,pim21,pim22,pim_dn,pim_up
2784  real(dp) :: pre11,pre12,pre21,pre22,pre_dn,pre_up
2785 !arrays
2786  real(dp) :: tsec(2)
2787 
2788 ! *************************************************************************
2789 
2790 !Real or complex inputs are coded
2791  DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
2792 
2793 !Real or complex outputs are coded
2794  DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot")
2795  DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
2796 
2797  DBG_CHECK(id >= 1,'Wrong id')
2798  DBG_CHECK(ip >= 1,'Wrong id')
2799 
2800  DBG_CHECK(multd >= 1,"wrong multd")
2801  DBG_CHECK(multp >= 1,"wrong multp")
2802 
2803  DBG_CHECK(nden-id-multd >=-1,'nden-id-multd')
2804  DBG_CHECK(npot-ip-multp >=-1,'npot-ip-multp')
2805 
2806  if(nspden==1)then
2807 
2808    if(cpldot==1 .or. cplex==1 )then
2809 
2810      do i2=1,multd
2811        do i1=1,multp
2812          ar=zero
2813 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar)
2814          do ir=1,cplex*nfft
2815            ar=ar + potarr(ir,1,ip+i1-1)*denarr(ir,1,id+i2-1)
2816          end do
2817          dot(1,i1,i2)=ar
2818        end do ! i1
2819      end do ! i2
2820 
2821    else  ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot
2822 
2823      do i2=1,multd
2824        do i1=1,multp
2825          ar=zero ; ai=zero
2826 !$OMP PARALLEL DO PRIVATE(ir,jr) SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai)
2827          do ir=1,nfft
2828            jr=2*ir
2829            ar=ar + potarr(jr-1,1,ip+i1-1)*denarr(jr-1,1,id+i2-1) &
2830 &           + potarr(jr  ,1,ip+i1-1)*denarr(jr  ,1,id+i2-1)
2831            ai=ai + potarr(jr-1,1,ip+i1-1)*denarr(jr  ,1,id+i2-1) &
2832 &           - potarr(jr  ,1,ip+i1-1)*denarr(jr-1,1,id+i2-1)
2833          end do
2834          dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai
2835        end do ! i1
2836      end do ! i2
2837 
2838    end if
2839 
2840  else if(nspden==2)then
2841 
2842    if(cpldot==1 .or. cplex==1 )then
2843 
2844      do i2=1,multd
2845        do i1=1,multp
2846          ar=zero
2847 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar)
2848          do ir=1,cplex*nfft
2849            ar=ar + potarr(ir,1,ip+i1-1)* denarr(ir,2,id+i2-1)               &       ! This is the spin up contribution
2850 &          + potarr(ir,2,ip+i1-1)*(denarr(ir,1,id+i2-1)-denarr(ir,2,id+i2-1)) ! This is the spin down contribution
2851          end do
2852          dot(1,i1,i2)=ar
2853        end do ! i1
2854      end do ! i2
2855 
2856    else ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot
2857 
2858      do i2=1,multd
2859        do i1=1,multp
2860          ar=zero ; ai=zero
2861 !$OMP PARALLEL DO PRIVATE(ir,jr,dre_up,dim_up,dre_dn,dim_dn,pre_up,pim_up,pre_dn,pim_dn) &
2862 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai)
2863          do ir=1,nfft
2864            jr=2*ir
2865 
2866            dre_up=denarr(jr-1,2,id+i2-1)
2867            dim_up=denarr(jr  ,2,id+i2-1)
2868            dre_dn=denarr(jr-1,1,id+i2-1)-dre_up
2869            dim_dn=denarr(jr  ,1,id+i2-1)-dim_up
2870 
2871            pre_up=potarr(jr-1,1,ip+i1-1)
2872            pim_up=potarr(jr  ,1,ip+i1-1)
2873            pre_dn=potarr(jr-1,2,ip+i1-1)
2874            pim_dn=potarr(jr  ,2,ip+i1-1)
2875 
2876            ar=ar + pre_up * dre_up &
2877 &           + pim_up * dim_up &
2878 &           + pre_dn * dre_dn &
2879 &           + pim_dn * dim_dn
2880            ai=ai + pre_up * dim_up &
2881 &           - pim_up * dre_up &
2882 &           + pre_dn * dim_dn &
2883 &           - pim_dn * dre_dn
2884 
2885          end do
2886          dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai
2887        end do ! i1
2888      end do ! i2
2889 
2890    end if
2891 
2892  else if(nspden==4)then
2893 !  \rho{\alpha,\beta} V^{\alpha,\beta} =
2894 !  rho*(V^{11}+V^{22})/2$
2895 !  + m_x Re(V^{12})- m_y Im{V^{12}}+ m_z(V^{11}-V^{22})/2
2896    if (cplex==1) then
2897      do i2=1,multd
2898        do i1=1,multp
2899          ar=zero
2900 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar)
2901          do ir=1,cplex*nfft
2902            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
2903 &          + potarr(ir,3,ip+i1-1)                                *denarr(ir,2,id+i2-1)& ! This is the m_x contrib
2904 &          - potarr(ir,4,ip+i1-1)                                *denarr(ir,3,id+i2-1)& ! This is the m_y contrib
2905 &          +(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
2906          end do
2907          dot(1,i1,i2)=ar
2908        end do ! i1
2909      end do ! i2
2910    else ! cplex=2
2911 !    Note concerning storage when cplex=2:
2912 !    V are stored as : v^11, v^22, V^12, i.V^21 (each are complex)
2913 !    N are stored as : n, m_x, m_y, mZ          (each are complex)
2914      if (cpldot==1) then
2915        do i2=1,multd
2916          do i1=1,multp
2917            ar=zero ; ai=zero
2918 !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre22,pim22,pre12,pim12) &
2919 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar)
2920            do ir=1,nfft
2921              jr=2*ir
2922              dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2))
2923              dim11=half*(denarr(jr  ,1,id+i2)+denarr(jr-1,4,id+i2))
2924              dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2))
2925              dim22=half*(denarr(jr  ,1,id+i2)-denarr(jr-1,4,id+i2))
2926              dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr  ,3,id+i2))
2927              dim12=half*(denarr(jr  ,2,id+i2)-denarr(jr-1,3,id+i2))
2928              dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr  ,3,id+i2))
2929              dim21=half*(denarr(jr  ,2,id+i2)+denarr(jr-1,3,id+i2))
2930              pre11= potarr(jr-1,1,ip+i1)
2931              pim11= potarr(jr  ,1,ip+i1)
2932              pre22= potarr(jr-1,2,ip+i1)
2933              pim22= potarr(jr  ,2,ip+i1)
2934              pre12= potarr(jr-1,3,ip+i1)
2935              pim12= potarr(jr  ,3,ip+i1)
2936              pre21= potarr(jr  ,4,ip+i1)
2937              pim21=-potarr(jr-1,4,ip+i1)
2938              ar=ar + pre11 * dre11 &
2939 &             + pim11 * dim11 &
2940 &             + pre22 * dre22 &
2941 &             + pim22 * dim22 &
2942 &             + pre12 * dre12 &
2943 &             + pim12 * dim12 &
2944 &             + pre21 * dre21 &
2945 &             + pim21 * dim21
2946            end do
2947            dot(1,i1,i2)=ar
2948          end do ! i1
2949        end do ! i2
2950      else !cpldot=2
2951        do i2=1,multd
2952          do i1=1,multp
2953            ar=zero ; ai=zero
2954 !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre12,pim12,pre22,pim22) &
2955 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai)
2956            do ir=1,nfft
2957              jr=2*ir
2958              dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2))
2959              dim11=half*(denarr(jr  ,1,id+i2)+denarr(jr-1,4,id+i2))
2960              dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2))
2961              dim22=half*(denarr(jr  ,1,id+i2)-denarr(jr-1,4,id+i2))
2962              dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr  ,3,id+i2))
2963              dim12=half*(denarr(jr  ,2,id+i2)-denarr(jr-1,3,id+i2))
2964              dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr  ,3,id+i2))
2965              dim21=half*(denarr(jr  ,2,id+i2)+denarr(jr-1,3,id+i2))
2966              pre11= potarr(jr-1,1,ip+i1)
2967              pim11= potarr(jr  ,1,ip+i1)
2968              pre22= potarr(jr-1,2,ip+i1)
2969              pim22= potarr(jr  ,2,ip+i1)
2970              pre12= potarr(jr-1,3,ip+i1)
2971              pim12= potarr(jr  ,3,ip+i1)
2972              pre21= potarr(jr  ,4,ip+i1)
2973              pim21=-potarr(jr-1,4,ip+i1)
2974              ar=ar + pre11 * dre11 &
2975 &             + pim11 * dim11 &
2976 &             + pre22 * dre22 &
2977 &             + pim22 * dim22 &
2978 &             + pre12 * dre12 &
2979 &             + pim12 * dim12 &
2980 &             + pre21 * dre21 &
2981 &             + pim21 * dim21
2982              ai=ai + pre11 * dim11 &
2983 &             - pim11 * dre11 &
2984 &             + pre22 * dim22 &
2985 &             - pim22 * dre22 &
2986 &             + pre12 * dim12 &
2987 &             - pim12 * dre12 &
2988 &             + pre21 * dim21 &
2989 &             - pim21 * dre21
2990            end do
2991            dot(1,i1,i2)=ar
2992            dot(2,i1,i2)=ai
2993          end do ! i1
2994        end do ! i2
2995      end if ! cpldot
2996    end if ! cplex
2997  end if ! nspden
2998 
2999  factor=ucvol/dble(nfftot)
3000  dot(:,:,:)=factor*dot(:,:,:)
3001 
3002 !XG030513 : MPIWF reduction (addition) on dot is needed here
3003  if (mpi_summarize) then
3004    call timab(48,1,tsec)
3005    call xmpi_sum(dot,mpicomm ,ierr)
3006    call timab(48,2,tsec)
3007  end if
3008 
3009  if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero
3010 
3011 end subroutine dotprodm_vn

ABINIT/findminscf [ Functions ]

[ Top ] [ Functions ]

NAME

 findminscf

FUNCTION

 Compute the minimum of a function whose value
 and derivative are known at two points, using different algorithms.
 Also deduce different quantities at this predicted
 point, and at the two other points

INPUTS

 choice=1,uses a linear interpolation of the derivatives
       =2,uses a quadratic interpolation based on the
        values of the function, and the second derivative at mid-point
 etotal_1=first value of the function
 etotal_2=second value of the function
 dedv_1=first value of the derivative
 dedv_2=second value of the derivative
 lambda_1=first value of the argument
 lambda_2=second value of the argument

OUTPUT

 dedv_predict=predicted value of the derivative (usually zero,
  except if choice=4, if it happens that a minimum cannot be located,
  and a trial step is taken)
 d2edv2_predict=predicted value of the second derivative (not if choice=4)
 d2edv2_1=first value of the second derivative (not if choice=4)
 d2edv2_2=second value of the second derivative (not if choice=4)
 etotal_predict=predicted value of the function
 lambda_predict=predicted value of the argument
 status= 0 if everything went normally ;
         1 if negative second derivative
         2 if some other problem

PARENTS

      scfcge

CHILDREN

      wrtout

SOURCE

2418 subroutine findminscf(choice,dedv_1,dedv_2,dedv_predict,&
2419 & d2edv2_1,d2edv2_2,d2edv2_predict,&
2420 & etotal_1,etotal_2,etotal_predict,&
2421 & lambda_1,lambda_2,lambda_predict,errid,errmess)
2422 
2423 
2424 !This section has been created automatically by the script Abilint (TD).
2425 !Do not modify the following lines by hand.
2426 #undef ABI_FUNC
2427 #define ABI_FUNC 'findminscf'
2428 !End of the abilint section
2429 
2430  implicit none
2431 
2432 !Arguments ------------------------------------
2433 !scalars
2434  integer,intent(in) :: choice
2435  integer,intent(out) :: errid
2436  character(len=500), intent(out) :: errmess
2437  real(dp),intent(in) :: dedv_1,dedv_2,etotal_1,etotal_2,lambda_1,lambda_2
2438  real(dp),intent(out) :: d2edv2_1,d2edv2_2,d2edv2_predict,dedv_predict
2439  real(dp),intent(out) :: etotal_predict,lambda_predict
2440 
2441 !Local variables-------------------------------
2442 !scalars
2443  real(dp) :: cc,d2edv2_mid,d_lambda,dedv_2bis
2444  real(dp) :: dedv_mid2,etotal_2bis
2445  character(len=500) :: message
2446 
2447 ! *************************************************************************
2448 
2449 !DEBUG
2450 !write(std_out,*)' findmin : enter'
2451 !write(std_out,*)' choice,lambda_1,lambda_2=',choice,lambda_1,lambda_2
2452 !ENDDEBUG
2453 
2454  errid = AB7_NO_ERROR
2455  d_lambda=lambda_1-lambda_2
2456 
2457  if(choice==1) then
2458 
2459 !  Use the derivative information to predict lambda
2460    d2edv2_mid=(dedv_1-dedv_2)/d_lambda
2461    lambda_predict=lambda_2-dedv_2/d2edv2_mid
2462    dedv_predict=dedv_2+(lambda_predict-lambda_2)*d2edv2_mid
2463    d2edv2_1=d2edv2_mid
2464    d2edv2_2=d2edv2_mid
2465    d2edv2_predict=d2edv2_mid
2466 !  also use the first energy to predict new energy
2467    etotal_predict=etotal_1+dedv_1*(lambda_predict-lambda_1)&
2468 &   +0.5_dp*d2edv2_1*(lambda_predict-lambda_1)**2
2469    etotal_2bis=etotal_1+dedv_1*(lambda_2-lambda_1)&
2470 &   +0.5_dp*d2edv2_1*(lambda_2-lambda_1)**2
2471 
2472    if(d2edv2_mid<0.0_dp)then
2473      errid = AB7_ERROR_MIXING_INTERNAL
2474      write(errmess,'(a,es18.10,a)')'The second derivative is negative, equal to ',d2edv2_mid,'.'
2475      MSG_WARNING(errmess)
2476    end if
2477 
2478  else if(choice==2) then
2479 
2480 !  Use energies and first derivative information
2481 !  etotal = aa + bb * lambda + cc * lambda**2
2482    dedv_mid2=(etotal_1-etotal_2)/d_lambda
2483    cc=(dedv_1-dedv_mid2)/d_lambda
2484    lambda_predict=lambda_1-0.5_dp*dedv_1/cc
2485    d2edv2_1=2*cc
2486    d2edv2_2=d2edv2_1
2487    d2edv2_predict=d2edv2_1
2488    if(d2edv2_predict<0.0_dp)then
2489      errid = AB7_ERROR_MIXING_INTERNAL
2490      write(errmess, '(a,es18.10,a,a,a)' )&
2491 &     'The second derivative is negative, equal to',d2edv2_predict,'.',ch10,&
2492 &     '=> Pivoting                     '
2493      MSG_WARNING(errmess)
2494      if(etotal_2 < etotal_1)then
2495        lambda_predict=lambda_2-0.5_dp*(lambda_1-lambda_2)
2496      else
2497        lambda_predict=lambda_1-0.5_dp*(lambda_2-lambda_1)
2498      end if
2499    end if
2500    dedv_predict=dedv_1+(lambda_predict-lambda_1)*d2edv2_1
2501    dedv_2bis=dedv_1+(lambda_2-lambda_1)*d2edv2_1
2502    etotal_predict=etotal_1+dedv_1*(lambda_predict-lambda_1)&
2503 &   +0.5_dp*d2edv2_1*(lambda_predict-lambda_1)**2
2504 
2505  end if
2506 
2507  write(message, '(a,es12.4,a,es18.10)' ) &
2508 & ' findmin : lambda_predict ',lambda_predict,' etotal_predict ',etotal_predict
2509  call wrtout(std_out,message,'COLL')
2510 
2511 end subroutine findminscf

ABINIT/m_ab7_mixing [ Modules ]

[ Top ] [ Modules ]

NAME

 m_ab7_mixing

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (XG, DC, 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 .

PARENTS

CHILDREN

SOURCE

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 
26 module m_ab7_mixing
27 
28  use defs_basis
29  use m_abicore
30  use m_errors
31  use m_linalg_interfaces
32  use m_xmpi
33 
34  use m_time,      only : timab
35  use m_io_tools,  only : open_file
36 
37  implicit none
38 
39  private

ABINIT/scfeig [ Functions ]

[ Top ] [ Functions ]

NAME

 scfeig

FUNCTION

 Compute the largest eigenvalue and eigenvector of the SCF cycle.
 A brute force algorithm is presently used.

INPUTS

  istep= number of the step in the SCF cycle
  nfft=(effective) number of FFT grid points (for this processor)
  nspden=number of spin-density components

OUTPUT

  (see side effects)

SIDE EFFECTS

  vtrial0(nfft,nspden)= contains vtrial at istep == 1
  vtrial(nfft,nspden)= at input, it is the trial potential that gave vresid .
       at output, it is an updated trial potential
  vrespc(nfft,nspden)=the input preconditioned residual potential
  work(nfft,nspden,2)=work space

PARENTS

      m_ab7_mixing

CHILDREN

      wrtout

SOURCE

1826 subroutine scfeig(istep,nfft,nspden,vrespc,vtrial,vtrial0,work,errid,errmess)
1827 
1828 
1829 !This section has been created automatically by the script Abilint (TD).
1830 !Do not modify the following lines by hand.
1831 #undef ABI_FUNC
1832 #define ABI_FUNC 'scfeig'
1833 !End of the abilint section
1834 
1835  implicit none
1836 
1837 !Arguments ------------------------------------
1838 !scalars
1839  integer,intent(in) :: istep,nfft,nspden
1840  integer,intent(out) :: errid
1841  character(len = 500), intent(out) :: errmess
1842 !arrays
1843  real(dp),intent(inout) :: vtrial0(nfft,nspden),work(nfft,nspden,2)
1844  real(dp),intent(inout) :: vrespc(nfft,nspden)
1845  real(dp), intent(inout) :: vtrial(nfft,nspden)
1846 
1847 !Local variables-------------------------------
1848 !scalars
1849  integer :: ifft,isp
1850  real(dp) :: eigen_scf,factor,fix_resid,resid_new,resid_old
1851  character(len=500) :: message
1852 
1853 ! *************************************************************************
1854 
1855  errid = AB7_NO_ERROR
1856 
1857  if(nspden==4)then
1858    errid = AB7_ERROR_MIXING_ARG
1859    write(errmess, *) ' scfeig : does not work yet for nspden=4'
1860    return
1861  end if
1862 
1863 !Set a fixed residual square for normalization of eigenvectors
1864  fix_resid=1.0d-4
1865 
1866 !A few initialisations for the first istep
1867  if(istep==1)then
1868 
1869    write(message, '(a,es12.4,a,a,a,a,a,a,a)' )&
1870 &   ' scfeig : fixed PC_residual square =',fix_resid,ch10,&
1871 &   '    Note that fixed resid should always be much larger',ch10,&
1872 &   '    than initial PC resid square, still sufficiently',ch10,&
1873 &   '    small to reduce anharmonic effects ',ch10
1874    call wrtout(std_out,message,'COLL')
1875 
1876 !  Compute the preconditioned residual
1877    resid_old=0.0_dp
1878    do isp=1,nspden
1879      do ifft=1,nfft
1880        resid_old=resid_old+vrespc(ifft,isp)**2
1881      end do
1882    end do
1883    write(message, '(a,es12.4)' )&
1884 &   ' scfeig : initial PC_residual square =',resid_old
1885    call wrtout(std_out,message,'COLL')
1886    if(resid_old>1.0d-8)then
1887      errid = AB7_ERROR_MIXING_ARG
1888      write(errmess,'(a,a,a,a,a,a,a,a,a,a)') ch10,&
1889 &     ' scfeig : ERROR -',ch10,&
1890 &     '  This value is not good enough to allow',ch10,&
1891 &     '  the computation of the eigenvectors of the SCF cycle.',ch10,&
1892 &     '  It should be better than 1.0d-8 .',ch10,&
1893 &     '  Action : improve the accuracy of your starting wavefunctions.'
1894      return
1895    end if
1896 
1897 !  Also transfer vtrial in vtrial_old
1898    vtrial0(:,:)=vtrial(:,:)
1899 
1900 !  In order to start the search for eigenvectors,
1901 !  use the tiny residual vector, renormalized
1902    factor=sqrt(fix_resid/resid_old)
1903    work(:,:,1)=vrespc(:,:)*factor
1904    vtrial(:,:)=vtrial0(:,:)+work(:,:,1)
1905 
1906 !  If istep is not equal to 1
1907  else if(istep>=2)then
1908 !
1909 !  Compute the corresponding operator expectation value
1910 !  And put the residual vector minus the difference
1911 !  between vtrial and vtrial_old
1912 !  (this is actually the action of the operator !) in vect(*,2)
1913    eigen_scf=0.0_dp
1914    do isp=1,nspden
1915      do ifft=1,nfft
1916        eigen_scf=eigen_scf+&
1917 &       work(ifft,isp,1) * vrespc(ifft,isp)
1918      end do
1919    end do
1920 
1921    do isp=1,nspden
1922      do ifft=1,nfft
1923        vrespc(ifft,isp)=vrespc(ifft,isp)&
1924 &       +vtrial(ifft,isp)-vtrial0(ifft,isp)
1925        work(ifft,isp,2)=vrespc(ifft,isp)
1926      end do
1927    end do
1928    eigen_scf=eigen_scf/fix_resid
1929    write(message, '(a,es12.4,a)' ) &
1930 &   ' scfeig : Operator expectation value ',eigen_scf,' (extremal eigenvalue * diemix)'
1931    call wrtout(std_out,message,'COLL')
1932    call wrtout(ab_out,message,'COLL')
1933 !
1934 !  Compute residual of vect(*,2)
1935    resid_new=zero
1936    do isp=1,min(nspden,2)
1937      do ifft=1,nfft
1938        resid_new=resid_new+ work(ifft,isp,2) ** 2
1939      end do
1940    end do
1941    if (nspden==4) then
1942      do ifft=1,nfft
1943        resid_new=resid_new+two*(work(ifft,3,2)**2+work(ifft,4,2)**2)
1944      end do
1945    end if
1946    factor=sqrt(fix_resid/resid_new)
1947    if(eigen_scf<zero) then
1948      factor=-factor ! the new vector MAY be oposite to the old one
1949 !    if(factor<-one) factor=-factor ! the new vector is not opposed to the old one
1950    end if
1951    write(message, '(a,es12.4)' ) &
1952 &   ' scfeig : Inverse of renormalization factor ',one/factor
1953    call wrtout(std_out,message,'COLL')
1954    call wrtout(ab_out,message,'COLL')
1955    write(message, '(a,es12.4)' ) &
1956 &   ' scfeig : Convergence criterion value (->0 at convergency) ',one/factor-eigen_scf-one
1957    call wrtout(std_out,message,'COLL')
1958    call wrtout(ab_out,message,'COLL')
1959 
1960    work(:,:,1)=work(:,:,2)*factor
1961    vtrial(:,:)=vtrial0(:,:)+work(:,:,1)
1962 !  End the different istep cases
1963  end if
1964 
1965 end subroutine scfeig

ABINIT/sqnormm_v [ Functions ]

[ Top ] [ Functions ]

NAME

 sqnormm_v

FUNCTION

 For a series of potentials,
 compute square of the norm (integral over FFT grid), to obtain
 a square residual-like quantity (so the sum of product of values
 is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume).
 Take into account the spin components of the density and potentials (nspden), and sum over them.
 Need the index of the first potential to be treated, in the provided array
 of potentials, and the number of potentials to be treated.
 Might be used to compute just one square of norm, in a big array, such as to avoid
 copying a potential from a big array to a temporary place.

INPUTS

  cplex=if 1, real space function on FFT grid is REAL, if 2, COMPLEX
  index=index of the first potential to be treated
  mpicomm=the mpi communicator used for the summation
  mpi_summarize=set it to .true. if parallelisation is done over FFT
  mult=number of potentials to be treated
  nfft= (effective) number of FFT grid points (for this processor)
  npot= third dimension of the potarr array
  nspden=number of spin-density components
  opt_storage: 0, if potential is stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn]
               1, if potential is stored as V, B_x, B_y, Bz  (B=magn. field)
  potarr(cplex*nfft,nspden,npot)=array of real space potentials on FFT grid

OUTPUT

  norm2(mult)= value of the square of the norm of the different potentials

PARENTS

      scfcge,scfopt

CHILDREN

      timab,xmpi_sum

SOURCE

3053 subroutine sqnormm_v(cplex,index,mpicomm, mpi_summarize,mult,nfft,norm2,npot,nspden,opt_storage,potarr)
3054 
3055 
3056 !This section has been created automatically by the script Abilint (TD).
3057 !Do not modify the following lines by hand.
3058 #undef ABI_FUNC
3059 #define ABI_FUNC 'sqnormm_v'
3060 !End of the abilint section
3061 
3062  implicit none
3063 
3064 !Arguments ------------------------------------
3065 !scalars
3066  integer,intent(in) :: cplex,index,mult,nfft,npot,nspden,opt_storage,mpicomm
3067  logical, intent(in) :: mpi_summarize
3068 !arrays
3069  real(dp),intent(in) :: potarr(cplex*nfft,nspden,npot)
3070  real(dp),intent(out) :: norm2(mult)
3071 
3072 !Local variables-------------------------------
3073 !scalars
3074  integer :: ierr,ifft,ii,ispden
3075  real(dp) :: ar
3076 !arrays
3077  real(dp) :: tsec(2)
3078 
3079 ! *************************************************************************
3080 
3081 !Real or complex inputs are coded
3082  DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
3083  DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
3084 
3085  DBG_CHECK(index>=1,"wrong index")
3086  DBG_CHECK(mult>=1,"wrong mult")
3087  DBG_CHECK(npot>=1,"wrong npot")
3088 
3089  DBG_CHECK(npot-index-mult>=-1,'npot-index-mult')
3090 
3091  do ii=1,mult
3092    ar=zero
3093    do ispden=1,min(nspden,2)
3094 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar)
3095      do ifft=1,cplex*nfft
3096        ar=ar + potarr(ifft,ispden,index+ii-1)**2
3097      end do
3098    end do
3099    norm2(ii)=ar
3100    if (nspden==4) then
3101      ar=zero
3102      do ispden=3,4
3103 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ii,index,ispden,nfft,potarr) REDUCTION(+:ar)
3104        do ifft=1,cplex*nfft
3105          ar=ar + potarr(ifft,ispden,index+ii-1)**2
3106        end do
3107      end do
3108      if (opt_storage==0) then
3109        if (cplex==1) then
3110          norm2(ii)=norm2(ii)+two*ar
3111        else
3112          norm2(ii)=norm2(ii)+ar
3113        end if
3114      else
3115        norm2(ii)=half*(norm2(ii)+ar)
3116      end if
3117    end if
3118  end do
3119 
3120 !XG030513 : MPIWF reduction (addition) on norm2 is needed here
3121  if (mpi_summarize) then
3122    call timab(48,1,tsec)
3123    call xmpi_sum(norm2,mpicomm ,ierr)
3124    call timab(48,2,tsec)
3125  end if
3126 
3127 end subroutine sqnormm_v

m_ab7_mixing/ab7_mixing_copy_current_step [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_copy_current_step

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      dfpt_newvtr,newrho,newvtr

CHILDREN

      nullify_

SOURCE

490 subroutine ab7_mixing_copy_current_step(mix, arr_resid, errid, errmess, &
491 &  arr_respc, arr_paw_resid, arr_paw_respc, arr_atm)
492 
493 !Arguments ------------------------------------
494 !scalars
495 
496 !This section has been created automatically by the script Abilint (TD).
497 !Do not modify the following lines by hand.
498 #undef ABI_FUNC
499 #define ABI_FUNC 'ab7_mixing_copy_current_step'
500 !End of the abilint section
501 
502  type(ab7_mixing_object), intent(inout) :: mix
503  real(dp), intent(in) :: arr_resid(mix%space * mix%nfft, mix%nspden)
504  integer, intent(out) :: errid
505  character(len = 500), intent(out) :: errmess
506  real(dp), intent(in), optional :: arr_respc(mix%space * mix%nfft, mix%nspden)
507  real(dp), intent(in), optional :: arr_paw_resid(mix%n_pawmix), arr_paw_respc(mix%n_pawmix)
508  real(dp), intent(in), optional :: arr_atm(3, mix%n_atom)
509 ! *************************************************************************
510 
511  if (.not. associated(mix%f_fftgr)) then
512     errid = AB7_ERROR_MIXING_ARG
513     write(errmess, '(a,a,a,a)' )ch10,&
514          & ' ab7_mixing_set_arr_current_step: ERROR -',ch10,&
515          & '  Working arrays not yet allocated.'
516     return
517  end if
518  errid = AB7_NO_ERROR
519 
520  mix%f_fftgr(:,:,mix%i_vresid(1)) = arr_resid(:,:)
521  if (present(arr_respc)) mix%f_fftgr(:,:,mix%i_vrespc(1)) = arr_respc(:,:)
522  if (present(arr_paw_resid)) mix%f_paw(:, mix%i_vresid(1)) = arr_paw_resid(:)
523  if (present(arr_paw_respc)) mix%f_paw(:, mix%i_vrespc(1)) = arr_paw_respc(:)
524  if (present(arr_atm)) mix%f_atm(:,:, mix%i_vresid(1)) = arr_atm(:,:)
525 
526 end subroutine ab7_mixing_copy_current_step

m_ab7_mixing/ab7_mixing_deallocate [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_deallocate

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      dfpt_scfcv,scfcv

CHILDREN

      nullify_

SOURCE

879 subroutine ab7_mixing_deallocate(mix)
880 
881 
882 !This section has been created automatically by the script Abilint (TD).
883 !Do not modify the following lines by hand.
884 #undef ABI_FUNC
885 #define ABI_FUNC 'ab7_mixing_deallocate'
886 !End of the abilint section
887 
888  implicit none
889 
890 !Arguments ------------------------------------
891 !scalars
892  type(ab7_mixing_object), intent(inout) :: mix
893 
894 !Local variables-------------------------------
895 !scalars
896  character(len = *), parameter :: subname = "ab7_mixing_deallocate"
897 ! *************************************************************************
898 
899  if (associated(mix%i_rhor)) then
900     ABI_DATATYPE_DEALLOCATE(mix%i_rhor)
901  end if
902  if (associated(mix%i_vtrial)) then
903     ABI_DATATYPE_DEALLOCATE(mix%i_vtrial)
904  end if
905  if (associated(mix%i_vresid)) then
906     ABI_DATATYPE_DEALLOCATE(mix%i_vresid)
907  end if
908  if (associated(mix%i_vrespc)) then
909     ABI_DATATYPE_DEALLOCATE(mix%i_vrespc)
910  end if
911  if (associated(mix%f_fftgr)) then
912     ABI_DEALLOCATE(mix%f_fftgr)
913  end if
914  if (associated(mix%f_paw)) then
915     ABI_DEALLOCATE(mix%f_paw)
916  end if
917  if (associated(mix%f_atm)) then
918     ABI_DEALLOCATE(mix%f_atm)
919  end if
920 
921  call nullify_(mix)
922 
923 end subroutine ab7_mixing_deallocate

m_ab7_mixing/ab7_mixing_eval [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_eval

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      dfpt_newvtr,newrho,newvtr

CHILDREN

      nullify_

SOURCE

716  subroutine ab7_mixing_eval(mix, arr, istep, nfftot, ucvol, &
717 & mpi_comm, mpi_summarize, errid, errmess, &
718 & reset, isecur, pawarr, pawopt, response, etotal, potden, &
719 & resnrm, comm_atom)
720 
721 
722 !This section has been created automatically by the script Abilint (TD).
723 !Do not modify the following lines by hand.
724 #undef ABI_FUNC
725 #define ABI_FUNC 'ab7_mixing_eval'
726 !End of the abilint section
727 
728  implicit none
729 
730 !Arguments ------------------------------------
731 !scalars
732  type(ab7_mixing_object), intent(inout) :: mix
733  integer, intent(in) :: istep, nfftot, mpi_comm
734  real(dp), intent(in) :: ucvol
735  real(dp), intent(inout) :: arr(mix%space * mix%nfft,mix%nspden)
736  logical, intent(in) :: mpi_summarize
737  integer, intent(out) :: errid
738  character(len = 500), intent(out) :: errmess
739  logical, intent(in), optional :: reset
740  integer, intent(in), optional :: isecur, comm_atom, pawopt, response
741  real(dp), intent(inout), optional :: pawarr(mix%n_pawmix)
742  real(dp), intent(in), optional :: etotal
743  real(dp), intent(in), optional :: potden(mix%space * mix%nfft,mix%nspden)
744  real(dp), intent(out), optional :: resnrm
745 
746 !Local variables-------------------------------
747 !scalars
748  integer :: moveAtm, dbl_nnsclo, initialized, isecur_
749  integer :: usepaw, pawoptmix_, response_
750  real(dp) :: resnrm_
751 ! *************************************************************************
752 
753  ! Argument checkings.
754  if (mix%iscf == AB7_MIXING_NONE) then
755     errid = AB7_ERROR_MIXING_ARG
756     write(errmess, '(a,a,a,a)' )ch10,&
757          & ' ab7_mixing_eval: ERROR -',ch10,&
758          & '  No method has been chosen.'
759     return
760  end if
761  if (mix%n_pawmix > 0 .and. .not. present(pawarr)) then
762     errid = AB7_ERROR_MIXING_ARG
763     write(errmess, '(a,a,a,a)' )ch10,&
764          & ' ab7_mixing_eval: ERROR -',ch10,&
765          & '  PAW is used, but no pawarr argument provided.'
766     return
767  end if
768  if (mix%n_atom > 0 .and. (.not. associated(mix%dtn_pc) .or. .not. associated(mix%xred))) then
769     errid = AB7_ERROR_MIXING_ARG
770     write(errmess, '(a,a,a,a)' )ch10,&
771          & ' ab7_mixing_eval: ERROR -',ch10,&
772          & '  Moving atoms is used, but no xred or dtn_pc attributes provided.'
773     return
774  end if
775  errid = AB7_NO_ERROR
776 
777  ! Miscellaneous
778  moveAtm = 0
779  if (mix%n_atom > 0) moveAtm = 1
780  initialized = 1
781  if (present(reset)) then
782     if (reset) initialized = 0
783  end if
784  isecur_ = 0
785  if (present(isecur)) isecur_ = isecur
786  usepaw = 0
787  if (mix%n_pawmix > 0) usepaw = 1
788  pawoptmix_ = 0
789  if (present(pawopt)) pawoptmix_ = pawopt
790  response_ = 0
791  if (present(response)) response_ = response
792 
793  ! Do the mixing.
794  resnrm_ = 0.d0
795  if (mix%iscf == AB7_MIXING_EIG) then
796     !  This routine compute the eigenvalues of the SCF operator
797     call scfeig(istep, mix%space * mix%nfft, mix%nspden, &
798          & mix%f_fftgr(:,:,mix%i_vrespc(1)), arr, &
799          & mix%f_fftgr(:,:,1), mix%f_fftgr(:,:,4:5), errid, errmess)
800  else if (mix%iscf == AB7_MIXING_SIMPLE .or. &
801       & mix%iscf == AB7_MIXING_ANDERSON .or. &
802       & mix%iscf == AB7_MIXING_ANDERSON_2 .or. &
803       & mix%iscf == AB7_MIXING_PULAY) then
804     if (present(comm_atom)) then
805       call scfopt(mix%space, mix%f_fftgr,mix%f_paw,mix%iscf,istep,&
806          & mix%i_vrespc,mix%i_vtrial, &
807          & mpi_comm,mpi_summarize,mix%nfft,mix%n_pawmix,mix%nspden, &
808          & mix%n_fftgr,mix%n_index,mix%kind,pawoptmix_,usepaw,pawarr, &
809          & resnrm_, arr, errid, errmess, comm_atom=comm_atom)
810     else
811       call scfopt(mix%space, mix%f_fftgr,mix%f_paw,mix%iscf,istep,&
812          & mix%i_vrespc,mix%i_vtrial, &
813          & mpi_comm,mpi_summarize,mix%nfft,mix%n_pawmix,mix%nspden, &
814          & mix%n_fftgr,mix%n_index,mix%kind,pawoptmix_,usepaw,pawarr, &
815          & resnrm_, arr, errid, errmess)
816     end if
817     !  Change atomic positions
818     if((istep==1 .or. mix%iscf==AB7_MIXING_SIMPLE) .and. mix%n_atom > 0)then
819        !    GAF: 2009-06-03
820        !    Apparently there are not reason
821        !    to restrict iscf=2 for ionmov=5
822        mix%xred(:,:) = mix%xred(:,:) + mix%dtn_pc(:,:)
823     end if
824  else if (mix%iscf == AB7_MIXING_CG_ENERGY .or.  mix%iscf == AB7_MIXING_CG_ENERGY_2) then
825     !  Optimize next vtrial using an algorithm based
826     !  on the conjugate gradient minimization of etotal
827     if (.not. present(etotal) .or. .not. present(potden)) then
828        errid = AB7_ERROR_MIXING_ARG
829        write(errmess, '(a,a,a,a)' )ch10,&
830             & ' ab7_mixing_eval: ERROR -',ch10,&
831             & '  Arguments etotal or potden are missing for CG on energy methods.'
832        return
833     end if
834     if (mix%n_atom == 0) then
835        ABI_ALLOCATE(mix%xred,(3,0))
836        ABI_ALLOCATE(mix%dtn_pc,(3,0))
837     end if
838     call scfcge(mix%space,dbl_nnsclo,mix%dtn_pc,etotal,mix%f_atm,&
839          & mix%f_fftgr,initialized,mix%iscf,isecur_,istep,&
840          & mix%i_rhor,mix%i_vresid,mix%i_vrespc,moveAtm,&
841          & mpi_comm,mpi_summarize,mix%n_atom,mix%nfft,nfftot,&
842          & mix%nspden,mix%n_fftgr,mix%n_index,mix%kind,&
843          & response_,potden,ucvol,arr,mix%xred, errid, errmess)
844     if (mix%n_atom == 0) then
845        ABI_DEALLOCATE(mix%xred)
846        ABI_DEALLOCATE(mix%dtn_pc)
847     end if
848     if (dbl_nnsclo == 1) errid = AB7_ERROR_MIXING_INC_NNSLOOP
849  end if
850 
851  if (present(resnrm)) resnrm = resnrm_
852 
853 end subroutine ab7_mixing_eval

m_ab7_mixing/ab7_mixing_eval_allocate [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_eval_allocate

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      dfpt_newvtr,newrho,newvtr

CHILDREN

      nullify_

SOURCE

552 subroutine ab7_mixing_eval_allocate(mix, istep)
553 
554 
555 !This section has been created automatically by the script Abilint (TD).
556 !Do not modify the following lines by hand.
557 #undef ABI_FUNC
558 #define ABI_FUNC 'ab7_mixing_eval_allocate'
559 !End of the abilint section
560 
561  implicit none
562 
563 !Arguments ------------------------------------
564 !scalars
565  type(ab7_mixing_object), intent(inout) :: mix
566  integer, intent(in), optional :: istep
567 
568 !Local variables-------------------------------
569 !scalars
570  integer :: istep_,temp_unit !, i_stat
571  real(dp) :: tsec(2)
572  character(len = *), parameter :: subname = "ab7_mixing_eval_allocate"
573  character(len=500) :: msg
574 
575 ! *************************************************************************
576 
577  istep_ = 1
578  if (present(istep)) istep_ = istep
579 
580  ! Allocate work array.
581  if (.not. associated(mix%f_fftgr)) then
582    !allocate(mix%f_fftgr(mix%space * mix%nfft,mix%nspden,mix%n_fftgr), stat = i_stat)
583    !call memocc_abi(i_stat, mix%f_fftgr, 'mix%f_fftgr', subname)
584    ABI_ALLOCATE(mix%f_fftgr,(mix%space * mix%nfft,mix%nspden,mix%n_fftgr))
585    mix%f_fftgr(:,:,:)=zero
586    if (mix%mffmem == 0 .and. istep_ > 1) then
587      call timab(83,1,tsec)
588      if (open_file(mix%diskCache,msg,newunit=temp_unit,form='unformatted',status='old') /= 0) then
589        MSG_ERROR(msg)
590      end if
591      rewind(temp_unit)
592      read(temp_unit) mix%f_fftgr
593      if (mix%n_pawmix == 0) close(unit=temp_unit)
594      call timab(83,2,tsec)
595    end if
596  end if
597  ! Allocate PAW work array.
598  if (.not. associated(mix%f_paw)) then
599     !allocate(mix%f_paw(mix%n_pawmix,mix%n_fftgr), stat = i_stat)
600     !call memocc_abi(i_stat, mix%f_paw, 'mix%f_paw', subname)
601     ABI_ALLOCATE(mix%f_paw,(mix%n_pawmix,mix%n_fftgr))
602     if (mix%n_pawmix > 0) then
603       mix%f_paw(:,:)=zero
604       if (mix%mffmem == 0 .and. istep_ > 1) then
605         read(temp_unit) mix%f_paw
606         close(unit=temp_unit)
607         call timab(83,2,tsec)
608       end if
609     end if
610  end if
611  ! Allocate atom work array.
612  if (.not. associated(mix%f_atm)) then
613     !allocate(mix%f_atm(3,mix%n_atom,mix%n_fftgr), stat = i_stat)
614     !call memocc_abi(i_stat, mix%f_atm, 'mix%f_atm', subname)
615     ABI_ALLOCATE(mix%f_atm,(3,mix%n_atom,mix%n_fftgr))
616  end if
617 
618  end subroutine ab7_mixing_eval_allocate

m_ab7_mixing/ab7_mixing_eval_deallocate [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_eval_deallocate

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      dfpt_newvtr,newrho,newvtr

CHILDREN

      nullify_

SOURCE

644  subroutine ab7_mixing_eval_deallocate(mix)
645 
646 
647 !This section has been created automatically by the script Abilint (TD).
648 !Do not modify the following lines by hand.
649 #undef ABI_FUNC
650 #define ABI_FUNC 'ab7_mixing_eval_deallocate'
651 !End of the abilint section
652 
653  implicit none
654 
655 !Arguments ------------------------------------
656 !scalars
657  type(ab7_mixing_object), intent(inout) :: mix
658 
659 !Local variables-------------------------------
660 !scalars
661  integer :: temp_unit !i_all, i_stat
662  real(dp) :: tsec(2)
663  character(len = *), parameter :: subname = "ab7_mixing_eval_deallocate"
664  character(len=500) :: msg
665 
666 ! *************************************************************************
667 
668  ! Save on disk and deallocate work array in case on disk cache only.
669  if (mix%mffmem == 0) then
670     call timab(83,1,tsec)
671     if (open_file(mix%diskCache,msg,newunit=temp_unit,form='unformatted',status='unknown') /= 0) then
672       MSG_ERROR(msg)
673     end if
674     rewind(temp_unit)
675     ! VALGRIND complains not all of f_fftgr_disk is initialized
676      write(temp_unit) mix%f_fftgr
677     if (mix%n_pawmix > 0) then
678       write(temp_unit) mix%f_paw
679     end if
680     close(unit=temp_unit)
681     call timab(83,2,tsec)
682     ABI_DEALLOCATE(mix%f_fftgr)
683     nullify(mix%f_fftgr)
684     if (associated(mix%f_paw)) then
685        ABI_DEALLOCATE(mix%f_paw)
686        nullify(mix%f_paw)
687     end if
688  end if
689 
690 end subroutine ab7_mixing_eval_deallocate

m_ab7_mixing/ab7_mixing_new [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_new

FUNCTION

INPUTS

OUTPUT

NOTES

PARENTS

      dfpt_scfcv,scfcv

CHILDREN

      nullify_

SOURCE

200 subroutine ab7_mixing_new(mix, iscf, kind, space, nfft, nspden, &
201 &  npawmix, errid, errmess, npulayit, useprec)
202 
203 
204 !This section has been created automatically by the script Abilint (TD).
205 !Do not modify the following lines by hand.
206 #undef ABI_FUNC
207 #define ABI_FUNC 'ab7_mixing_new'
208 !End of the abilint section
209 
210  implicit none
211 
212 !Arguments ------------------------------------
213 !scalars
214  type(ab7_mixing_object), intent(out) :: mix
215  integer, intent(in) :: iscf, kind, space, nfft, nspden, npawmix
216  integer, intent(out) :: errid
217  character(len = 500), intent(out) :: errmess
218  integer, intent(in), optional :: npulayit
219  logical, intent(in), optional :: useprec
220 
221 !Local variables-------------------------------
222 !scalars
223  integer :: ii !, i_stat
224  character(len = *), parameter :: subname = "ab7_mixing_new"
225 ! *************************************************************************
226 
227  ! Set default values.
228  call init_(mix)
229 
230  ! Argument checkings.
231  if (kind /= AB7_MIXING_POTENTIAL .and. kind /= AB7_MIXING_DENSITY) then
232     errid = AB7_ERROR_MIXING_ARG
233     write(errmess, '(a,a,a,a)' )ch10,&
234          & ' ab7_mixing_set_arrays: ERROR -',ch10,&
235          & '  Mixing must be done on density or potential only.'
236     return
237  end if
238  if (space /= AB7_MIXING_REAL_SPACE .and. space /= AB7_MIXING_FOURRIER_SPACE) then
239     errid = AB7_ERROR_MIXING_ARG
240     write(errmess, '(a,a,a,a)' )ch10,&
241          & ' ab7_mixing_set_arrays: ERROR -',ch10,&
242          & '  Mixing must be done in real or Fourrier space only.'
243     return
244  end if
245  if (iscf /= AB7_MIXING_EIG .and. iscf /= AB7_MIXING_SIMPLE .and. &
246       & iscf /= AB7_MIXING_ANDERSON .and. &
247       & iscf /= AB7_MIXING_ANDERSON_2 .and. &
248       & iscf /= AB7_MIXING_CG_ENERGY .and. &
249       & iscf /= AB7_MIXING_PULAY .and. &
250       & iscf /= AB7_MIXING_CG_ENERGY_2) then
251     errid = AB7_ERROR_MIXING_ARG
252     write(errmess, "(A,I0,A)") "Unknown mixing scheme (", iscf, ")."
253     return
254  end if
255  errid = AB7_NO_ERROR
256 
257  ! Mandatory arguments.
258  mix%iscf     = iscf
259  mix%kind     = kind
260  mix%space    = space
261  mix%nfft     = nfft
262  mix%nspden   = nspden
263  mix%n_pawmix = npawmix
264 
265  ! Optional arguments.
266  if (present(useprec)) mix%useprec = useprec
267 
268  ! Set-up internal dimensions.
269  !These arrays are needed only in the self-consistent case
270  if (iscf == AB7_MIXING_EIG) then
271     !    For iscf==1, five additional vectors are needed
272     !    The index 1 is attributed to the old trial potential,
273     !    The new residual potential, and the new
274     !    preconditioned residual potential receive now a temporary index
275     !    The indices number 4 and 5 are attributed to work vectors.
276     mix%n_fftgr=5 ; mix%n_index=1
277  else if(iscf == AB7_MIXING_SIMPLE) then
278     !    For iscf==2, three additional vectors are needed.
279     !    The index number 1 is attributed to the old trial vector
280     !    The new residual potential, and the new preconditioned
281     !    residual potential, receive now a temporary index.
282     mix%n_fftgr=3 ; mix%n_index=1
283     if (.not. mix%useprec) mix%n_fftgr = 2
284  else if(iscf == AB7_MIXING_ANDERSON) then
285     !    For iscf==3 , four additional vectors are needed.
286     !    The index number 1 is attributed to the old trial vector
287     !    The new residual potential, and the new and old preconditioned
288     !    residual potential, receive now a temporary index.
289     mix%n_fftgr=4 ; mix%n_index=2
290     if (.not. mix%useprec) mix%n_fftgr = 3
291  else if (iscf == AB7_MIXING_ANDERSON_2) then
292     !    For iscf==4 , six additional vectors are needed.
293     !    The indices number 1 and 2 are attributed to two old trial vectors
294     !    The new residual potential, and the new and two old preconditioned
295     !    residual potentials, receive now a temporary index.
296     mix%n_fftgr=6 ; mix%n_index=3
297     if (.not. mix%useprec) mix%n_fftgr = 5
298  else if(iscf == AB7_MIXING_CG_ENERGY .or. iscf == AB7_MIXING_CG_ENERGY_2) then
299     !    For iscf==5 or 6, ten additional vectors are needed
300     !    The index number 1 is attributed to the old trial vector
301     !    The index number 6 is attributed to the search vector
302     !    Other indices are attributed now. Altogether ten vectors
303     mix%n_fftgr=10 ; mix%n_index=3
304  else if(iscf == AB7_MIXING_PULAY) then
305     !    For iscf==7, lot of additional vectors are needed
306     !    The index number 1 is attributed to the old trial vector
307     !    The index number 2 is attributed to the old residual
308     !    The indices number 2 and 3 are attributed to two old precond. residuals
309     !    Other indices are attributed now.
310     if (present(npulayit)) mix%n_pulayit = npulayit
311     mix%n_fftgr=2+2*mix%n_pulayit ; mix%n_index=1+mix%n_pulayit
312     if (.not. mix%useprec) mix%n_fftgr = 1+2*mix%n_pulayit
313  end if ! iscf cases
314 
315  ! Allocate new arrays.
316  !allocate(mix%i_rhor(mix%n_index), stat = i_stat)
317  !call memocc_abi(i_stat, mix%i_rhor, 'mix%i_rhor', subname)
318  ABI_DATATYPE_ALLOCATE(mix%i_rhor,(mix%n_index))
319  mix%i_rhor(:)=0
320  !allocate(mix%i_vtrial(mix%n_index), stat = i_stat)
321  !call memocc_abi(i_stat, mix%i_vtrial, 'mix%i_vtrial', subname)
322  ABI_DATATYPE_ALLOCATE(mix%i_vtrial,(mix%n_index))
323  mix%i_vtrial(:)=0
324  !allocate(mix%i_vresid(mix%n_index), stat = i_stat)
325  !call memocc_abi(i_stat, mix%i_vresid, 'mix%i_vresid', subname)
326  ABI_DATATYPE_ALLOCATE(mix%i_vresid,(mix%n_index))
327  mix%i_vresid(:)=0
328  !allocate(mix%i_vrespc(mix%n_index), stat = i_stat)
329  !call memocc_abi(i_stat, mix%i_vrespc, 'mix%i_vrespc', subname)
330  ABI_DATATYPE_ALLOCATE(mix%i_vrespc,(mix%n_index))
331  mix%i_vrespc(:)=0
332 
333  ! Setup initial values.
334  if (iscf == AB7_MIXING_EIG) then
335     mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2 ; mix%i_vrespc(1)=3
336  else if(iscf == AB7_MIXING_SIMPLE) then
337     mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2 ; mix%i_vrespc(1)=3
338     if (.not. mix%useprec) mix%i_vrespc(1)=2
339  else if(iscf == AB7_MIXING_ANDERSON) then
340     mix%i_vtrial(1)=1 ; mix%i_vresid(1)=2
341     if (mix%useprec) then
342        mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=4
343     else
344        mix%i_vrespc(1)=2 ; mix%i_vrespc(2)=3
345     end if
346  else if (iscf == AB7_MIXING_ANDERSON_2) then
347     mix%i_vtrial(1)=1 ; mix%i_vtrial(2)=2
348     mix%i_vresid(1)=3
349     if (mix%useprec) then
350        mix%i_vrespc(1)=4 ; mix%i_vrespc(2)=5 ; mix%i_vrespc(3)=6
351     else
352        mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=4 ; mix%i_vrespc(3)=5
353     end if
354  else if(iscf == AB7_MIXING_CG_ENERGY .or. iscf == AB7_MIXING_CG_ENERGY_2) then
355     mix%n_fftgr=10 ; mix%n_index=3
356     mix%i_vtrial(1)=1
357     mix%i_vresid(1)=2 ; mix%i_vresid(2)=4 ; mix%i_vresid(3)=7
358     mix%i_vrespc(1)=3 ; mix%i_vrespc(2)=5 ; mix%i_vrespc(3)=8
359     mix%i_rhor(2)=9 ; mix%i_rhor(3)=10
360  else if(iscf == AB7_MIXING_PULAY) then
361     do ii=1,mix%n_pulayit
362        mix%i_vtrial(ii)=2*ii-1 ; mix%i_vrespc(ii)=2*ii
363     end do
364     mix%i_vrespc(mix%n_pulayit+1)=2*mix%n_pulayit+1
365     mix%i_vresid(1)=2*mix%n_pulayit+2
366     if (.not. mix%useprec) mix%i_vresid(1)=2
367  end if ! iscf cases
368 
369 end subroutine ab7_mixing_new

m_ab7_mixing/ab7_mixing_use_disk_cache [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_use_disk_cache

FUNCTION

INPUTS

OUTPUT

NOTES

  Obsolete?

PARENTS

      dfpt_scfcv,scfcv

CHILDREN

      nullify_

SOURCE

393 subroutine ab7_mixing_use_disk_cache(mix, fnametmp_fft)
394 
395 
396 !This section has been created automatically by the script Abilint (TD).
397 !Do not modify the following lines by hand.
398 #undef ABI_FUNC
399 #define ABI_FUNC 'ab7_mixing_use_disk_cache'
400 !End of the abilint section
401 
402  implicit none
403 
404 !Arguments ------------------------------------
405 !scalars
406  type(ab7_mixing_object), intent(inout) :: mix
407  character(len = *), intent(in) :: fnametmp_fft
408 ! *************************************************************************
409 
410  if (len(trim(fnametmp_fft)) > 0) then
411     mix%mffmem = 0
412     write(mix%diskCache, "(A)") fnametmp_fft
413  else
414     mix%mffmem = 1
415  end if
416 
417 end subroutine ab7_mixing_use_disk_cache

m_ab7_mixing/ab7_mixing_use_moving_atoms [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  ab7_mixing_use_moving_atoms

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      newrho,newvtr

CHILDREN

      nullify_

SOURCE

443 subroutine ab7_mixing_use_moving_atoms(mix, natom, xred, dtn_pc)
444 
445 !Arguments ------------------------------------
446 !scalars
447 
448 !This section has been created automatically by the script Abilint (TD).
449 !Do not modify the following lines by hand.
450 #undef ABI_FUNC
451 #define ABI_FUNC 'ab7_mixing_use_moving_atoms'
452 !End of the abilint section
453 
454  type(ab7_mixing_object), intent(inout) :: mix
455  integer, intent(in) :: natom
456  real(dp), intent(in), target :: dtn_pc(3, natom)
457  real(dp), intent(in), target :: xred(3, natom)
458 
459 ! *************************************************************************
460 
461  mix%n_atom = natom
462  mix%dtn_pc => dtn_pc
463  mix%xred => xred
464 
465 end subroutine ab7_mixing_use_moving_atoms

m_ab7_mixing/init_ [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  init_

FUNCTION

  Initialize the object

PARENTS

      m_ab7_mixing

CHILDREN

      nullify_

SOURCE

106 subroutine init_(mix)
107 
108 
109 !This section has been created automatically by the script Abilint (TD).
110 !Do not modify the following lines by hand.
111 #undef ABI_FUNC
112 #define ABI_FUNC 'init_'
113 !End of the abilint section
114 
115  implicit none
116 
117 !Arguments ------------------------------------
118 !scalars
119  type(ab7_mixing_object), intent(out) :: mix
120 ! *************************************************************************
121 
122  ! Default values.
123  mix%iscf      = AB7_MIXING_NONE
124  mix%mffmem    = 1
125  mix%n_index   = 0
126  mix%n_fftgr   = 0
127  mix%n_pulayit = 7
128  mix%n_pawmix  = 0
129  mix%n_atom    = 0
130  mix%useprec   = .true.
131 
132  call nullify_(mix)
133 
134 end subroutine init_

m_ab7_mixing/nullify [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

  nullify_

FUNCTION

  Nullify the pointers

PARENTS

      m_ab7_mixing

CHILDREN

      nullify_

SOURCE

152 subroutine nullify_(mix)
153 
154 
155 !This section has been created automatically by the script Abilint (TD).
156 !Do not modify the following lines by hand.
157 #undef ABI_FUNC
158 #define ABI_FUNC 'nullify_'
159 !End of the abilint section
160 
161  implicit none
162 
163 !Arguments ------------------------------------
164 !scalars
165  type(ab7_mixing_object), intent(inout) :: mix
166 ! *************************************************************************
167 
168  ! Nullify internal pointers.
169  nullify(mix%i_rhor)
170  nullify(mix%i_vtrial)
171  nullify(mix%i_vresid)
172  nullify(mix%i_vrespc)
173  nullify(mix%f_fftgr)
174  nullify(mix%f_atm)
175  nullify(mix%f_paw)
176 
177 end subroutine nullify_

m_ab7_mixing/scfcge [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

 scfcge

FUNCTION

 Compute the next vtrial of the SCF cycle.
 Uses a conjugate gradient minimization of the total energy
 Can move only the trial potential (if moved_atm_inside==0), or
 move the trial atomic positions as well (if moved_atm_inside==1).

INPUTS

  cplex= if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
  dtn_pc(3,natom)=preconditioned change of atomic position, in reduced
    coordinates. Will be quickly transferred to f_atm(:,:,i_vrespc(1))
  etotal=the actual total energy
  initialized= if 0, the initialization of the gstate run is not yet finished
  iscf =5 => SCF cycle, CG based on estimation of energy gradient
       =6 => SCF cycle, CG based on true minimization of the energy
  isecur=level of security of the computation
  istep= number of the step in the SCF cycle
  moved_atm_inside: if==1, the atoms are allowed to move.
  mpicomm=the mpi communicator used for the summation
  mpi_summarize=set it to .true. if parallelisation is done over FFT
  natom=number of atoms
  nfft=(effective) number of FFT grid points (for this processor)
  nfftot=total number of FFT grid points
  nspden=number of spin-density components
  n_fftgr=third dimension of the array f_fftgr
  n_index=dimension for indices of potential/density (see i_vresid, ivrespc, i_rhor...)
  opt_denpot= 0 vtrial (and also f_fftgr) really contains the trial potential
              1 vtrial (and also f_fftgr) actually contains the trial density
  response= if 0, GS calculation, if 1, RF calculation, intrinsically harmonic !
  rhor(cplex*nfft,nspden)=actual density
  ucvol=unit cell volume in bohr**3

OUTPUT

 dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.

SIDE EFFECTS

 Input/Output:
  vtrial(cplex*nfft,nspden)= at input, it is the trial potential that gave
       the input residual of the potential and Hellman-Feynman forces
                       at output, it is the new trial potential .
  xred(3,natom)=(needed if moved_atm_inside==1)
      reduced dimensionless atomic coordinates
      at input, those that generated the input residual of the potential
      and Hellman-Feynman forces, at output, these are the new ones.
  f_fftgr(cplex*nfft,nspden,n_fftgr)=different functions defined on the fft grid :
   The input vtrial is transferred, at output, in f_fftgr(:,:,1).
   The input f_fftgr(:,:,i_vresid(1)) contains the last residual.
     the value of i_vresid(1) is transferred to i_vresid(2) at output.
   The input f_fftgr(:,:,i_vresid(2)) contains the old residual.
     the value of i_vresid(2) is transferred to i_vresid(3) at output.
   The input f_fftgr(:,:,i_vresid(3)) contains the previous last residual.
   For the preconditioned potential residual, the same logic as for the
     the potential residual is used, with i_vrespc replacing i_vresid.
   The input rhor is transferred, at output, in f_fft(:,:,i_rhor(2)).
   The old density is input in f_fft(:,:,i_rhor(2)), and the value of
      i_rhor(2) is transferred to i_rhor(3) before the end of the routine.
   The input/output search vector is stored in f_fftgr(:,:,6)
  f_atm(3,natom,n_fftgr)=different functions defined for each atom :
   The input xred is transferred, at output, in f_atm(:,:,1).
   The input f_atm(:,:,i_vresid(1)) contains minus the HF forces.
     the value of i_vresid(1) is transferred to i_vresid(2) at output.
   The input f_atm(:,:,i_vresid(2)) contains minus the old HF forces.
     the value of i_vresid(2) is transferred to i_vresid(3) at output.
   The input f_atm(:,:,i_vresid(3)) contains minus the previous old HF forces.
   For the preconditioned change of atomic positions, the same logic as for the
     the potential residual is used, with i_vrespc replacing i_vresid.
   The input/output search vector is stored in f_atm(:,:,6)
  i_rhor(2:3)=index of the density (past and previous past) in the array f_fftgr
  i_vresid(3)=index of the residual potentials (present, past and previous
   past) in the array f_fftgr; also similar index for minus Hellman-Feynman
   forces in the array f_atm .
  i_vrespc(3)=index of the preconditioned residual potentials
                  (present, past and previous past) in the array f_fftgr ;
   also similar index for the preconditioned change of atomic position (dtn_pc).

TODO

 This routine is much too difficult to read ! Should be rewritten ...
 Maybe make separate subroutines for line search and CG step ?!

PARENTS

      m_ab7_mixing

CHILDREN

      aprxdr,findminscf,sqnormm_v,wrtout

SOURCE

1017 subroutine scfcge(cplex,dbl_nnsclo,dtn_pc,etotal,f_atm,&
1018 & f_fftgr,initialized,iscf,isecur,istep,&
1019 & i_rhor,i_vresid,i_vrespc,moved_atm_inside,mpicomm,mpi_summarize,&
1020 & natom,nfft,nfftot,nspden,n_fftgr,n_index,opt_denpot,response,rhor,ucvol,vtrial,xred,errid,errmess)
1021 
1022 
1023 !This section has been created automatically by the script Abilint (TD).
1024 !Do not modify the following lines by hand.
1025 #undef ABI_FUNC
1026 #define ABI_FUNC 'scfcge'
1027 !End of the abilint section
1028 
1029  implicit none
1030 
1031 !Arguments ------------------------------------
1032 !scalars
1033  integer,intent(in) :: cplex,initialized,iscf,isecur,istep,moved_atm_inside,mpicomm
1034  integer,intent(in) :: n_fftgr,n_index,natom,nfft,nfftot,nspden,opt_denpot,response
1035  integer,intent(out) :: dbl_nnsclo, errid
1036  character(len = 500), intent(out) :: errmess
1037  logical, intent(in) :: mpi_summarize
1038  real(dp),intent(in) :: etotal,ucvol
1039 !arrays
1040  integer,intent(inout) :: i_rhor(n_index),i_vresid(n_index),i_vrespc(n_index)
1041  real(dp),intent(in) :: dtn_pc(3,natom),rhor(cplex*nfft,nspden)
1042  real(dp),intent(inout) :: f_atm(3,natom,n_fftgr)
1043  real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr)
1044  real(dp),intent(inout) :: vtrial(cplex*nfft,nspden),xred(3,natom)
1045 
1046 !Local variables-------------------------------
1047 !mlinmin gives the maximum number of steps in the line minimization
1048 !   after which the algorithm is restarted (with a decrease of the
1049 !   adaptative trial step length). This number should not be large,
1050 !   since if the potential landscape is harmonic, the number of
1051 !   search steps should be small. If it is large, we are not in the
1052 !   harmonic region, and the CG algorithm will not be really useful,
1053 !   so one can just restart the algorithm ...
1054 !scalars
1055  integer,parameter :: mlinmin=5
1056  integer,save :: end_linmin,iline_cge,ilinear,ilinmin,isecur_eff,nlinear
1057  integer,save :: number_of_restart,status
1058  integer :: choice,iatom,idir,ifft,iline_cge_input,ilinmin_input,isp
1059  integer :: testcg,tmp,errid_
1060  real(dp),save :: d2edv2_old2,d_lambda_old2,dedv_old2,etotal_old
1061  real(dp),save :: etotal_previous,lambda_adapt,lambda_new,lambda_old,resid_old
1062  real(dp) :: d2e11,d2e12,d2e22,d2edv2_new,d2edv2_old
1063  real(dp) :: d2edv2_predict,d_lambda,de1,de2,dedv_mix
1064  real(dp) :: dedv_new,dedv_old,dedv_predict,determ,etotal_input
1065  real(dp) :: etotal_predict,gamma,lambda_input,lambda_predict2
1066  real(dp) :: lambda_predict=1.0_dp,ratio,reduction
1067  real(dp) :: resid_input,temp
1068  character(len=500) :: message
1069 !arrays
1070  real(dp) :: resid_new(1)
1071  real(dp), allocatable :: tmp_fft1(:,:)
1072 
1073 ! *************************************************************************
1074 
1075  errid = AB7_NO_ERROR
1076  dbl_nnsclo = 0
1077 
1078 !reduction gives the level of reduction of the error in
1079 !the line minimization to be reached for the minimization to be
1080 !considered successfull
1081  reduction=0.1_dp
1082 
1083 !nlinear increases with the number of times the 2D minimization succeded
1084 !to reach the true minimum directly. It is a measure of the
1085 !degree of parabolicity of the problem, and is used to
1086 !skip some steps by performing extrapolation.
1087  if(istep==1)then
1088 
1089 !  Skipping some steps is sometimes unsecure, so it is possible
1090 !  to make nlinear start at a negative value - if isecur is positive
1091    isecur_eff=isecur
1092    nlinear=min(-isecur_eff,0)
1093    ilinear=0
1094 
1095 !  Response function calculation are intrinsically harmonic, so one
1096 !  can shift isecur (by -2), and start with a positive nlinear
1097    if(response==1)then
1098      isecur_eff=isecur-2
1099      nlinear=-isecur_eff
1100      ilinear=nlinear
1101    end if
1102 
1103    iline_cge=0
1104    ilinmin=0
1105  end if
1106 
1107 !Compute actual residual resid_new (residual of f_fftgr(:,:,i_vrespc(1))
1108  call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr)
1109 
1110 !Save input residual and ilinmin for final printing
1111  resid_input=resid_new(1)
1112  etotal_input=etotal
1113  ilinmin_input=ilinmin
1114  iline_cge_input=iline_cge
1115 !Transfer dtn_pc in f_atm
1116  if(moved_atm_inside==1)then
1117    f_atm(:,:,i_vrespc(1))=dtn_pc(:,:)
1118  end if
1119 
1120 !=======================================================================
1121 !Now the routine is decomposed in three mutually exclusive parts :
1122 !if(istep==1)then initialize the algorithm
1123 !else if(ilinmin>0)then perform the line minimisation
1124 !else if(ilinmin==0)then determine the new search direction (CG step)
1125 !=======================================================================
1126 
1127 
1128 !--------------------------------------
1129 !Here initialize the algorithm
1130  if(istep==1)then
1131 
1132 !  At the beginning of each gstate run, lambda_adapt is forced to have the
1133 !  same value, that is 1.0_dp. In the other cases when istep=1 (at different
1134 !  broyden steps, for example), the previously obtained
1135 !  adaptive value is kept.
1136    if(initialized==0)lambda_adapt=1.0_dp
1137    lambda_old=0.0_dp
1138    lambda_input=0.0_dp
1139    number_of_restart=0
1140    lambda_new=lambda_adapt
1141 
1142    f_fftgr(:,:,1)=vtrial(:,:)
1143    f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1144 
1145 !  This copy must be written in F77, because of stack problems on the DECs
1146    do isp=1,nspden
1147      do ifft=1,cplex*nfft
1148        f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1))
1149      end do
1150    end do
1151    vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6)
1152    if(moved_atm_inside==1)then
1153      f_atm(:,:,1)=xred(:,:)
1154      f_atm(:,:,i_rhor(2))=xred(:,:)
1155 !    There shouldn t be problems with the stack size for this small array.
1156      f_atm(:,:,6)=f_atm(:,:,i_vrespc(1))
1157      xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6)
1158    end if
1159    tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
1160    tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1161    ilinmin=1
1162    resid_old=resid_new(1)
1163    etotal_old=etotal
1164 
1165    status=0
1166 
1167 !  --------------------------------------
1168 
1169 !  Here performs the line minimisation
1170  else if(ilinmin>0)then
1171 
1172    lambda_input=lambda_new
1173 
1174 !  The choice with the Brent algorithm has been abandoned in version 1.6.m
1175 
1176 !  Compute the approximate energy derivatives dedv_new and dedv_old,
1177 !  from vresid and vresid_old
1178    choice=2
1179    call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,&
1180 &   f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,&
1181 &   natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred)
1182    d_lambda=lambda_new-lambda_old
1183    dedv_old=dedv_old/d_lambda
1184    dedv_new=dedv_new/d_lambda
1185 
1186 !  DEBUG
1187 !  write(std_out,'(a,4es12.4,i3)' )' scfcge:lold,lnew,dold,dnew,status',  &
1188 !  &  lambda_old,lambda_new,dedv_old,dedv_new,status
1189 !  ENDDEBUG
1190 
1191    if(status==0 .or. status==3)then
1192 !
1193 !    Then, compute a predicted point along the line
1194 !    The value of choice determines the minimization algorithm
1195 !    choice=1 uses the two values of the derivative of the energy
1196 !    choice=2 uses the two values of the energy, and and estimate of the
1197 !    second derivative at the mid-point.
1198 
1199      choice=1
1200      if(iscf==6)choice=2
1201      call findminscf(choice,dedv_new,dedv_old,dedv_predict,&
1202 &     d2edv2_new,d2edv2_old,d2edv2_predict,&
1203 &     etotal,etotal_old,etotal_predict,&
1204 &     lambda_new,lambda_old,lambda_predict,errid_,message)
1205      if (errid_ /= AB7_NO_ERROR) then
1206        call wrtout(std_out,message,'COLL')
1207      end if
1208 
1209 !    Suppress the next line for debugging  (there is another such line)
1210      status=0
1211 
1212 !    DEBUG
1213 !    Keep this debugging feature : it gives access to the investigation of lines
1214 !    in a different approach
1215 !    if(response==1 .and. istep>8)then
1216 !    lambda_predict=1.2d-2
1217 !    if(istep>=15)lambda_predict=lambda_predict-0.002
1218 !    if(istep>=14)stop
1219 !    status=3
1220 !    end if
1221 !    ENDDEBUG
1222 
1223    else
1224      if(status/=-1)then
1225        status=-1
1226        lambda_predict=-2.5_dp
1227      else
1228        lambda_predict=lambda_predict+0.1_dp
1229      end if
1230    end if
1231 
1232 !  If the predicted point is very close to the most recent
1233 !  computed point, while this is the first trial on this line,
1234 !  then we are in the linear regime :
1235 !  nlinear is increased by one unit. For the time being, do this even when
1236 !  moved_atm_inside==1 (the code still works when it is done, but it
1237 !  seems to be a bit unstable). The maximal value of nlinear is 1, except
1238 !  when isecur_eff is a negative number, less than -1.
1239    if( abs(lambda_predict-lambda_new)/&
1240 &   (abs(lambda_predict)+abs(lambda_new)) < 0.01 .and. ilinmin==1  ) then
1241 !    if(moved_atm_inside==0 .and. nlinear<max(1,-isecur_eff) )nlinear=nlinear+1
1242      if(nlinear<max(1,-isecur_eff))nlinear=nlinear+1
1243      ilinear=nlinear
1244    end if
1245 
1246 !  If the predicted point is close to the most recent computed point,
1247 !  or the previous one, set on the flag of end of line minization
1248    end_linmin=0
1249    if(abs(lambda_new-lambda_predict)*2.0_dp&
1250 &   /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1
1251    if(abs(lambda_old-lambda_predict)*2.0_dp&
1252 &   /(abs(lambda_predict)+abs(lambda_new)) <reduction) end_linmin=1
1253 
1254    if(status/=0)end_linmin=0
1255 
1256 !  Save the closest old lambda, if needed,
1257 !  also examine the reduction of the interval, and eventual stop
1258 !  the present line minimisation, because of convergence (end_linmin=1)
1259 !  Also treat the case in which the predicted value of lambda is negative,
1260 !  or definitely too small in which case the algorithm has to be restarted
1261 !  (not a very good solution, though ...)
1262 !  Finally also treat the case where insufficiently converged
1263 !  density at lambda=0.0_dp happens, which screws up the line minimisation.
1264 
1265 !  Here restart the algorithm with the best vtrial.
1266 !  Also make reduction in lambda_adapt
1267 !  DEBUG
1268 !  write(std_out,*)' scfcge : status=',status
1269 !  ENDDEBUG
1270    if( end_linmin==0 .and. status==0 .and.                               &
1271 &   (  (lambda_predict<0.005_dp*lambda_adapt .and. iscf==5)     .or.  &
1272 &   (abs(lambda_predict)<0.005_dp*lambda_adapt .and. iscf==6).or.  &
1273 &   ilinmin==mlinmin                                      )     )then
1274      if(number_of_restart>12)then
1275        errid = AB7_ERROR_MIXING_CONVERGENCE
1276        write(errmess,'(a,a,i0,a,a,a,a,a)')&
1277 &       'Potential-based CG line minimization not',' converged after ',number_of_restart,' restarts. ',ch10,&
1278 &       'Action : read the eventual warnings about lack of convergence.',ch10,&
1279 &       'Some might be relevant. Otherwise, raise nband. Returning'
1280        MSG_WARNING(errmess)
1281        return
1282      end if
1283 !    Make reduction in lambda_adapt (kind of steepest descent...)
1284      write(message,'(a,a,a)')&
1285 &     'Potential-based CG line minimization has trouble to converge.',ch10,&
1286 &     'The algorithm is restarted with more secure parameters.'
1287      MSG_WARNING(message)
1288      number_of_restart=number_of_restart+1
1289 !    At the second restart, double the number of non-self consistent loops.
1290      if(number_of_restart>=2)dbl_nnsclo=1
1291      lambda_adapt=lambda_adapt*0.7_dp
1292      lambda_new=lambda_adapt
1293 !    If the last energy is better than the old one, transfer the data.
1294 !    Otherwise, no transfer must occur (very simple to code...)
1295      if(etotal<etotal_old .or. abs(lambda_old)<1.0d-8)then
1296        f_fftgr(:,:,1)=vtrial(:,:)
1297        f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1298        do isp=1,nspden
1299          do ifft=1,cplex*nfft
1300            f_fftgr(ifft,isp,6)=f_fftgr(ifft,isp,i_vrespc(1))
1301          end do
1302        end do
1303        if(moved_atm_inside==1)then
1304          f_atm(:,:,1)=xred(:,:)
1305          f_atm(:,:,i_rhor(2))=xred(:,:)
1306          f_atm(:,:,6)=f_atm(:,:,i_vrespc(1))
1307        end if
1308        tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
1309        tmp=i_vresid(2) ; i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1310        resid_old=resid_new(1)
1311        etotal_old=etotal
1312      end if
1313      lambda_old=0.0_dp
1314      ilinmin=1
1315 !    Putting the flag to -1 avoids the usual actions taken with end_linmin=1
1316      end_linmin=-1
1317 !    Also put ilinear and nlinear to 0
1318      ilinear=0
1319      nlinear=0
1320 
1321 !    Here lambda_new is the closest to lambda_predict,
1322 !    or lambda_old is still 0.0_dp, while the energy shows that the minimum
1323 !    is away from 0.0_dp (insufficiently converged density at lambda=0.0_dp).
1324    else if( abs(lambda_new-lambda_predict)<abs(lambda_old-lambda_predict) &
1325 &     .or.                                                           &
1326 &     ( abs(lambda_old)<1.0d-6 .and.                               &
1327 &     ilinmin>1              .and.                               &
1328 &     etotal>etotal_previous         )                           &
1329 &     )then
1330      f_fftgr(:,:,1)=vtrial(:,:)
1331      tmp=i_rhor(3) ; i_rhor(3)=i_rhor(2) ; i_rhor(2)=tmp
1332      f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1333      tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2)
1334      i_vrespc(2)=i_vrespc(1); i_vrespc(1)=tmp;
1335      tmp=i_vresid(3); i_vresid(3)=i_vresid(2)
1336      i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1337      if(moved_atm_inside==1)then
1338        f_atm(:,:,1)=xred(:,:)
1339        f_atm(:,:,i_rhor(2))=xred(:,:)
1340      end if
1341      d_lambda_old2=lambda_old-lambda_new
1342      lambda_old=lambda_new
1343      etotal_old=etotal
1344      resid_old=resid_new(1)
1345      d2edv2_old2=d2edv2_new
1346      dedv_old=dedv_new
1347      dedv_old2=dedv_new
1348 !    if(abs(lambda_new-lambda_predict)*2.0_dp&
1349 !    &    /abs(lambda_new+lambda_predict)        <reduction) end_linmin=1
1350 
1351 !    Here lambda_old is the closest to lambda_predict (except for avoiding
1352 !    lambda_old==0.0_dp)
1353    else
1354      tmp=i_vresid(3) ; i_vresid(3)=i_vresid(1) ; i_vresid(1)=tmp
1355      f_fftgr(:,:,i_rhor(3))=rhor(:,:)
1356      if(moved_atm_inside==1) f_atm(:,:,i_rhor(3))=xred(:,:)
1357      tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(1) ; i_vrespc(1)=tmp
1358      d_lambda_old2=lambda_new-lambda_old
1359      etotal_previous=etotal
1360      d2edv2_old2=d2edv2_old
1361      dedv_old2=dedv_old
1362 !    if(abs(lambda_old-lambda_predict)*2.0_dp&
1363 !    &    /abs(lambda_old+lambda_predict)        <reduction) end_linmin=1
1364    end if
1365 
1366 !  If the interval has not yet been sufficiently reduced,
1367 !  continue the search
1368    if(end_linmin==0)then
1369      lambda_new=lambda_predict
1370 
1371 !    DEBUG
1372 !    write(std_out,'(a,2es16.6)' )&
1373 !    &   ' scfcge : continue search, lambda_old,lambda_new=',lambda_old,lambda_new
1374 !    write(std_out,'(a,2es16.6)' )&
1375 !    &   ' scfcge : f_fftgr(3:4,1,1)=',f_fftgr(3:4,1,1)
1376 !    write(std_out,'(a,2es16.6)' )&
1377 !    &   ' scfcge : f_fftgr(3:4,1,6)=',f_fftgr(3:4,1,6)
1378 !    ENDDEBUG
1379 
1380      vtrial(:,:)=f_fftgr(:,:,1)+(lambda_new-lambda_old)*f_fftgr(:,:,6)
1381      if(moved_atm_inside==1)then
1382        xred(:,:)=f_atm(:,:,1)+(lambda_new-lambda_old)*f_atm(:,:,6)
1383      end if
1384 
1385      ilinmin=ilinmin+1
1386 !
1387 !    Here generates a starting point for next line search
1388    else
1389      iline_cge=iline_cge+1
1390      if(end_linmin==1)ilinmin=0
1391      lambda_old=0.0_dp
1392 
1393 !    In order to generate the new step, take into account previous
1394 !    optimal lambdas (including those of previous ion moves),
1395 !    and the selected new one, if it is positive.
1396 !    However, wait iline_cge>1 to select new ones.
1397 !    lambda_adapt has been initialized at 1.0_dp
1398      if(iline_cge>1 .and. lambda_new>0.0_dp )then
1399 !      Actually compute a geometric mean
1400        lambda_adapt= ( lambda_adapt**(dble(iline_cge-1)) * abs(lambda_new)) &
1401 &       **(1.0_dp/dble(iline_cge))
1402 !      In order to recover the previous algorithm, it is enough
1403 !      to decomment the next line
1404 !      lambda_adapt=1.0_dp
1405      end if
1406      lambda_new=lambda_adapt
1407 
1408      vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2))
1409      if(moved_atm_inside==1)then
1410        xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2))
1411      end if
1412 
1413 !    End choice between continue line minim and determine new direction
1414    end if
1415 
1416 !
1417 !  -------------------------------
1418 
1419 !  Here perform the CG step
1420 
1421  else if(ilinmin==0)then
1422 
1423 !  Compute the approximate energy derivatives dedv_mix,dedv_new,dedv_old
1424    choice=3
1425    call aprxdr(cplex,choice,dedv_mix,dedv_new,dedv_old,&
1426 &   f_atm,f_fftgr,i_rhor(2),i_vresid,moved_atm_inside,mpicomm,mpi_summarize,&
1427 &   natom,nfft,nfftot,nspden,n_fftgr,rhor,ucvol,xred)
1428 
1429    dedv_mix=dedv_mix/lambda_new
1430    dedv_new=dedv_new/lambda_new
1431    dedv_old=dedv_old/lambda_new
1432 
1433 !  DEBUG
1434 !  write(message, '(a,3es12.4)' )' scfcge: lambda_adapt',&
1435 !  &     lambda_adapt
1436 !  call wrtout(std_out,message,'COLL')
1437 
1438 !  write(message, '(a,3es12.4)' )' scfcge: dedv_old,dedv_new,dedv_mix',&
1439 !  &     dedv_old,dedv_new,dedv_mix
1440 !  call wrtout(std_out,message,'COLL')
1441 !  ENDDEBUG
1442 
1443 !  Then, compute a predicted point, either along the line,
1444 !  or in a 2D plane
1445    testcg=1
1446    if(testcg==0)then
1447 !    This part corresponds to steepest descent,
1448 !    in which the line minimisation can be done
1449 !    using different algorithms, varying with the value of choice
1450      choice=1
1451      if(iscf==6)choice=2
1452      call findminscf(choice,dedv_new,dedv_old,dedv_predict,&
1453 &     d2edv2_new,d2edv2_old,d2edv2_predict,&
1454 &     etotal,etotal_old,etotal_predict,&
1455 &     lambda_new,lambda_old,lambda_predict,errid_,message)
1456      if (errid_ /= AB7_NO_ERROR) then
1457        call wrtout(std_out,message,'COLL')
1458      end if
1459      lambda_predict2=0.0_dp
1460 !    Suppress the next line for debugging (there is another such line)
1461      status=0
1462    else
1463 !    This part corresponds to conjugate gradient
1464 !    A 2D minimisation is performed
1465 !    oldest direction is labelled 2
1466 !    newest direction is labelled 1
1467      de1=dedv_old ;  de2=dedv_old2
1468      d2e11=(dedv_new-dedv_old)/lambda_new
1469      d2e22=d2edv2_old2
1470      d2e12=(dedv_mix-dedv_old)/d_lambda_old2
1471 !    The system to be solved is
1472 !    0 = de1 + lambda1 d2e11 + lambda2 d2d12
1473 !    0 = de2 + lambda1 d2e12 + lambda2 d2d22
1474      determ=d2e11*d2e22-d2e12*d2e12
1475      lambda_predict=-(de1*d2e22-de2*d2e12)/determ
1476      lambda_predict2=(de1*d2e12-de2*d2e11)/determ
1477      d2edv2_new=d2e11 ;  d2edv2_old=d2e11
1478    end if
1479 
1480 !  DEBUG
1481 !  write(message, '(a,5es11.3)' )' scfcge: de1,de2,d2e11,d2e22,d2e12',&
1482 !  &               de1,de2,d2e11,d2e22,d2e12
1483 !  call wrtout(std_out,message,'COLL')
1484 !  write(std_out,'(a,2es12.4)' )' scfcge: la_predict,la_predict2',&
1485 !  &               lambda_predict,lambda_predict2
1486 !  -----
1487 !  write(std_out,*)'residues ',
1488 !  !$       de1+lambda_predict*d2e11+lambda_predict2*d2e12,
1489 !  !$       de2+lambda_predict*d2e12+lambda_predict2*d2e22
1490 !  if(.true.)stop
1491 !  ENDDEBUG
1492 !
1493 
1494 !  Determine the region of the 2D search space
1495 !  in which the predicted point is located,
1496 !  or use linear indicator to decide interpolation
1497 !  and advance to next 2D search.
1498    end_linmin=0
1499    write(message, '(a,2i3)' )' nlinear, ilinear',nlinear,ilinear
1500    call wrtout(std_out,message,'COLL')
1501    if(lambda_predict<0.0_dp)then
1502 !    Something is going wrong. Just take a reasonable step
1503 !    along the steepest descent direction (Region III).
1504 !    Actually, Region I and region III are treated in the same way later.
1505 !    In effect, this corresponds to restart the algorithm
1506      end_linmin=3
1507 !    Also put ilinear and nlinear to 0
1508      ilinear=0
1509      nlinear=0
1510 !    Decrease the adaptive step to predict next direction
1511      lambda_adapt=lambda_adapt*0.7_dp
1512    else if(ilinear>=1) then
1513 !    Region IV : will do an interpolation
1514      end_linmin=4
1515      ilinear=ilinear-1
1516    else if(abs(lambda_predict2)>reduction          .or.&
1517 &     lambda_predict<0.5_dp                .or.&
1518 &     lambda_predict>2.5_dp                .or.&
1519 &     lambda_predict-abs(lambda_predict2)/reduction <0.0_dp  ) then
1520 !    Region II : lambda_predict is not too good, and not too bad.
1521      end_linmin=2
1522    else if (abs(1.0_dp-lambda_predict)<reduction)then
1523 !    Region I, the out-of-line point is OK.
1524      end_linmin=1
1525    else
1526 !    If everything fails, then region II.
1527      end_linmin=2
1528    end if
1529 
1530 !  DEBUG
1531 !  write(message, '(a,2es12.4,i2)' )&
1532 !  &     ' scfcge : la_predict, la_predict2, region',&
1533 !  &       lambda_predict,lambda_predict2,end_linmin
1534 !  call wrtout(std_out,message,'COLL')
1535 !  ENDDEBUG
1536 
1537 !  Treat region I, in the same way as region III
1538    if(end_linmin==1 .or. end_linmin==3)then
1539 
1540 !    In region I, the line search is
1541 !    along vtrial-vtrial_old.
1542 !    The closest point is the new point
1543 !    thus to be transfered in the "old" locations
1544 
1545      do isp=1,nspden
1546        do ifft=1,cplex*nfft
1547          f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new
1548        end do
1549      end do
1550      f_fftgr(:,:,1)=vtrial(:,:)
1551      f_fftgr(:,:,i_rhor(2))=rhor(:,:)
1552      if(moved_atm_inside==1)then
1553        f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new
1554        f_atm(:,:,1)=xred(:,:)
1555        f_atm(:,:,i_rhor(2))=xred(:,:)
1556      end if
1557      tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
1558      tmp=i_vresid(3) ; i_vresid(3)=i_vresid(2)
1559      i_vresid(2)=i_vresid(1) ; i_vresid(1)=tmp
1560      d_lambda_old2=-lambda_new
1561      lambda_old=lambda_new
1562      etotal_old=etotal
1563      resid_old=resid_new(1)
1564      d2edv2_old=d2edv2_new
1565      dedv_old=dedv_new
1566 
1567 !    Region I or III : one is close of the 2D minimum,
1568 !    or lambda_predict was negative (indicate a problem of convergence)
1569 !    Compute next trial potential along the
1570 !    PC residual and not along this search direction.
1571      ilinmin=0
1572 !    Question : isn t it here that one should prevent region I to called
1573 !    itself more than 1 time ???
1574 !    Here the small difference between region I and region III
1575      if(end_linmin==3)ilinmin=1
1576      lambda_old=0.0_dp
1577      lambda_new=lambda_adapt
1578 
1579      vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2))
1580      if(moved_atm_inside==1)then
1581        xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2))
1582      end if
1583 !    The new vtrial has been generated
1584 
1585    else
1586 
1587 !    Here region II or IV
1588      ilinmin=1
1589      if (lambda_predict==0._dp) then
1590        gamma=zero
1591      else
1592        gamma=lambda_predict2/lambda_predict
1593      end if
1594 !    Compute new search direction and trial potential
1595      write(message,*)' compute new search direction '
1596      call wrtout(std_out,message,'COLL')
1597      do isp=1,nspden
1598        do ifft=1,cplex*nfft
1599          f_fftgr(ifft,isp,6)=(vtrial(ifft,isp)-f_fftgr(ifft,isp,1))/lambda_new+ &
1600 &         gamma*f_fftgr(ifft,isp,6)
1601        end do
1602      end do
1603      vtrial(:,:)=f_fftgr(:,:,1)+ lambda_predict*f_fftgr(:,:,6)
1604      if(moved_atm_inside==1)then
1605        f_atm(:,:,6)=(xred(:,:)-f_atm(:,:,1))/lambda_new+ gamma*f_atm(:,:,6)
1606        xred(:,:)=f_atm(:,:,1)+ lambda_predict*f_atm(:,:,6)
1607      end if
1608 
1609 !    If end_linmin==2, then this vtrial is the good one
1610 
1611      if(end_linmin==2)then
1612 
1613        lambda_old=0.0_dp
1614        lambda_new=lambda_predict
1615 
1616      else if(end_linmin==4)then
1617 
1618 !      predict the result of the computation at the trial potential
1619 !      defined in the end_linmin==2 case
1620        gamma=lambda_predict2/d_lambda_old2
1621        ratio=lambda_predict/lambda_new
1622 
1623 !      Take care of vtrial
1624        f_fftgr(:,:,1)=vtrial(:,:)
1625 
1626        ABI_ALLOCATE(tmp_fft1,(cplex*nfft,nspden))
1627 !      Take care of vresid
1628        tmp_fft1(:,:)=f_fftgr(:,:,i_vresid(2))
1629        f_fftgr(:,:,i_vresid(2))=tmp_fft1(:,:)&
1630 &       +ratio*(f_fftgr(:,:,i_vresid(1))-tmp_fft1(:,:))&
1631 &       +gamma*(f_fftgr(:,:,i_vresid(3))-tmp_fft1(:,:))
1632        f_fftgr(:,:,i_vresid(3))=tmp_fft1(:,:)
1633 
1634 !      Take care of rhor
1635        tmp_fft1(:,:)=f_fftgr(:,:,i_rhor(2))
1636        f_fftgr(:,:,i_rhor(2))=tmp_fft1(:,:)&
1637 &       +ratio*(rhor(:,:)-tmp_fft1(:,:))&
1638 &       +gamma*(f_fftgr(:,:,i_rhor(3))-tmp_fft1(:,:))
1639        f_fftgr(:,:,i_rhor(3))=tmp_fft1(:,:)
1640 
1641 !      Take care of vrespc
1642        tmp_fft1(:,:)=f_fftgr(:,:,i_vrespc(2))
1643        f_fftgr(:,:,i_vrespc(2))=tmp_fft1(:,:)&
1644 &       +ratio*(f_fftgr(:,:,i_vrespc(1))-tmp_fft1(:,:))&
1645 &       +gamma*(f_fftgr(:,:,i_vrespc(3))-tmp_fft1(:,:))
1646        f_fftgr(:,:,i_vrespc(3))=tmp_fft1(:,:)
1647        ABI_DEALLOCATE(tmp_fft1)
1648 
1649        if(moved_atm_inside==1)then
1650          do idir=1,3
1651            do iatom=1,natom
1652 
1653 !            Take care of xred
1654              f_atm(idir,iatom,1)=xred(idir,iatom)
1655 
1656 !            Take care of -HF forces
1657              temp=f_atm(idir,iatom,i_vresid(2))
1658              f_atm(idir,iatom,i_vresid(2))=f_atm(idir,iatom,i_vresid(2))&
1659 &             +ratio*(f_atm(idir,iatom,i_vresid(1))-f_atm(idir,iatom,i_vresid(2)))&
1660 &             +gamma*(f_atm(idir,iatom,i_vresid(3))-f_atm(idir,iatom,i_vresid(2)))
1661              f_atm(idir,iatom,i_vresid(3))=temp
1662 
1663 !            Take care of old xreds
1664              temp=f_atm(idir,iatom,i_rhor(2))
1665              f_atm(idir,iatom,i_rhor(2))=f_atm(idir,iatom,i_rhor(2))&
1666 &             +ratio*(   xred(idir,iatom)          -f_atm(idir,iatom,i_rhor(2)))&
1667 &             +gamma*(f_atm(idir,iatom,i_rhor(3))-f_atm(idir,iatom,i_rhor(2)))
1668              f_atm(idir,iatom,i_rhor(3))=temp
1669 
1670 !            Take care of preconditioned changes of atomic positions
1671              temp=f_atm(idir,iatom,i_vrespc(2))
1672              f_atm(idir,iatom,i_vrespc(2))=f_atm(idir,iatom,i_vrespc(2))&
1673 &             +ratio*(f_atm(idir,iatom,i_vrespc(1))-f_atm(idir,iatom,i_vrespc(2)))&
1674 &             +gamma*(f_atm(idir,iatom,i_vrespc(3))-f_atm(idir,iatom,i_vrespc(2)))
1675              f_atm(idir,iatom,i_vrespc(3))=temp
1676 
1677            end do
1678          end do
1679        end if
1680 
1681 !      Since we are at the 2D minimum, the derivative is supposed
1682 !      to vanish. Note that dedv_old should not change, by contrast.
1683        dedv_old2=0.0_dp
1684        d_lambda_old2=-lambda_predict
1685        d2edv2_old2=-dedv_old/lambda_predict
1686        lambda_old=lambda_predict
1687        ilinmin=0
1688 
1689 !      So, jump to the next line
1690        iline_cge=iline_cge+1
1691        write(message,*)' energy CG update : after 2D interpolation,'
1692        call wrtout(std_out,message,'COLL')
1693        write(message,*)'    computation in the next plane '
1694        call wrtout(std_out,message,'COLL')
1695        write(message,*)
1696        call wrtout(std_out,message,'COLL')
1697        lambda_old=0.0_dp
1698        lambda_new=lambda_adapt
1699 
1700        vtrial(:,:)=f_fftgr(:,:,1)+lambda_new*f_fftgr(:,:,i_vrespc(2))
1701        if(moved_atm_inside==1)then
1702          xred(:,:)=f_atm(:,:,1)+lambda_new*f_atm(:,:,i_vrespc(2))
1703        end if
1704 
1705 !      The new trial potential is now generated
1706 
1707 !      End the specific treatment of region IV
1708      end if
1709 !
1710 !    End the choice between treatment of region I, II, or IV
1711    end if
1712 
1713 !  End of choice between initialisation or more developed parts of the CG algorithm
1714  else
1715    errid = AB7_ERROR_MIXING_ARG
1716    errmess = 'scfcge : BUG You should not be here ! '
1717    return
1718  end if
1719 
1720 !--------------------------------------
1721 
1722 !Write information : it will be easy to read by typing  grep scfcge logfile
1723 
1724  if(istep==1)then
1725    write(message,'(a,a,a)') ' scfcge:',ch10,' scfcge:istep-iline_cge-ilinmin lambda      etot             resid '
1726    call wrtout(std_out,message,'COLL')
1727  end if
1728 
1729  if(ilinmin_input/=0 .or. istep==1)then
1730 !  Usual line minimisation step
1731 
1732    if(iline_cge_input<10)then
1733      write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' )&
1734 &     ' scfcge: actual  ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input
1735    else
1736      write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' )&
1737 &     ' scfcge: actual  ',istep,'-',iline_cge_input,'-',ilinmin_input,lambda_input,etotal_input,resid_input
1738    end if
1739    call wrtout(std_out,message,'COLL')
1740 
1741    if( (end_linmin==1.or.end_linmin==-1) .and. istep/=1 )then
1742 
1743      if(end_linmin==1)then
1744        write(message, '(a,es13.4,a,i2,a,a)' )&
1745 &       ' scfcge: predict         ',lambda_predict,&
1746 &       ' suff. close => next line, ilinear=',ilinear,ch10,&
1747 &       ' scfcge:'
1748      else if(end_linmin==-1)then
1749        write(message, '(a,es13.4,a,a,a)' )&
1750 &       ' scfcge: predict         ',lambda_predict,&
1751 &       ' restart the algorithm ',ch10,&
1752 &       ' scfcge:'
1753      end if
1754      call wrtout(std_out,message,'COLL')
1755 
1756      if(iline_cge_input<9)then
1757        write(message, '(a,i4,a,i1,a,i1,es13.4,es20.12,es12.4)' ) &
1758 &       ' scfcge: start   ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old
1759      else
1760        write(message, '(a,i3,a,i2,a,i1,es13.4,es20.12,es12.4)' ) &
1761 &       ' scfcge: start   ',istep,'-',iline_cge,'-',0,0.0,etotal_old,resid_old
1762      end if
1763      call wrtout(std_out,message,'COLL')
1764 
1765    else if(istep/=1) then
1766      write(message, '(a,es13.4,a)' )&
1767 &     ' scfcge: predict         ',lambda_predict,&
1768 &     ' not close enough => continue minim.'
1769      call wrtout(std_out,message,'COLL')
1770    end if
1771 
1772  else
1773 !  CG prediction
1774    if(iline_cge_input<10)then
1775      write(message, '(a,i4,a,i1,a,es11.4,es20.12,es12.4,a,i1)' )&
1776 &     ' scfcge: actual  ',istep,'-',iline_cge_input,'-off',&
1777 &     lambda_adapt,etotal_input,resid_input,', end=',end_linmin
1778    else
1779      write(message, '(a,i3,a,i2,a,es11.4,es20.12,es12.4,a,i1)' )&
1780 &     ' scfcge: actual  ',istep,'-',iline_cge_input,'-off',&
1781 &     lambda_adapt,etotal_input,resid_input,', end=',end_linmin
1782    end if
1783    call wrtout(std_out,message,'COLL')
1784 
1785    if(end_linmin==4)then
1786      write(message, '(a)' ) ' scfcge:'
1787      call wrtout(std_out,message,'COLL')
1788    end if
1789 
1790  end if
1791 
1792 end subroutine scfcge

m_ab7_mixing/scfopt [ Functions ]

[ Top ] [ m_ab7_mixing ] [ Functions ]

NAME

 scfopt

FUNCTION

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

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

2027 subroutine scfopt(cplex,f_fftgr,f_paw,iscf,istep,i_vrespc,i_vtrial,&
2028 & mpicomm,mpi_summarize,nfft,npawmix,nspden,n_fftgr,&
2029 & n_index,opt_denpot,pawoptmix,usepaw,vpaw,vresid,vtrial,errid,errmess, &
2030 & comm_atom) ! optional
2031 
2032 
2033 !This section has been created automatically by the script Abilint (TD).
2034 !Do not modify the following lines by hand.
2035 #undef ABI_FUNC
2036 #define ABI_FUNC 'scfopt'
2037 !End of the abilint section
2038 
2039  implicit none
2040 
2041 !Arguments ------------------------------------
2042 !scalars
2043  integer,intent(in) :: cplex,iscf,istep,n_fftgr,n_index,nfft
2044  integer,intent(in) :: npawmix,nspden,opt_denpot,pawoptmix,usepaw,mpicomm
2045  integer, intent(in),optional :: comm_atom
2046  integer,intent(out) :: errid
2047  character(len = 500), intent(out) :: errmess
2048  logical, intent(in) :: mpi_summarize
2049  real(dp), intent(out) :: vresid
2050 !arrays
2051  integer,intent(inout) :: i_vrespc(n_index),i_vtrial(n_index)
2052  real(dp),intent(inout) :: f_fftgr(cplex*nfft,nspden,n_fftgr)
2053  real(dp),intent(inout) :: f_paw(npawmix,n_fftgr*usepaw),vpaw(npawmix*usepaw)
2054  real(dp),intent(inout) :: vtrial(cplex*nfft,nspden)
2055 !Local variables-------------------------------
2056 !scalars
2057  integer,parameter :: npulaymax=50
2058  integer :: i_vstore,ierr,ifft,ii,index,isp,jj,comm_atom_,niter,npulay,tmp
2059  real(dp),save :: prod_resid_old,resid_old,resid_old2
2060  real(dp) :: aa1,aa2,bb,cc1,cc2,current,det,lambda,lambda2,resid_best
2061  character(len=500) :: message
2062 !arrays
2063  integer,allocatable :: ipiv(:)
2064  real(dp),save :: amat(npulaymax+1,npulaymax+1)
2065  real(dp) :: mpibuff(2),prod_resid(1),prod_resid2(1),resid_new(1)
2066  real(dp),allocatable :: alpha(:),amatinv(:,:),amat_paw(:),rwork(:)
2067 
2068 ! *************************************************************************
2069 
2070 !DEBUG
2071 !write(std_out,*)' scfopt : enter ; istep,iscf ',istep,iscf
2072 !ENDDEBUG
2073 
2074  errid = AB7_NO_ERROR
2075 
2076  comm_atom_=xmpi_comm_self; if(present(comm_atom)) comm_atom_=comm_atom
2077 
2078  i_vstore=i_vtrial(1)
2079  if (iscf==4) i_vstore=i_vtrial(2)
2080  if (iscf==7) then
2081    if (modulo(n_fftgr, 2) == 0 ) then
2082      npulay=(n_fftgr-2)/2
2083    else
2084      npulay=(n_fftgr-1)/2
2085    end if
2086    i_vstore=i_vtrial(npulay)
2087  else
2088    npulay=0
2089  end if
2090 
2091 !Compute the new residual resid_new, from f_fftgr/f_paw(:,:,i_vrespc(1))
2092  call sqnormm_v(cplex,i_vrespc(1),mpicomm,mpi_summarize,1,nfft,resid_new,n_fftgr,nspden,opt_denpot,f_fftgr)
2093  if (usepaw==1.and.pawoptmix==1) then
2094    do index=1,npawmix
2095      resid_new(1)=resid_new(1)+f_paw(index,i_vrespc(1))**2
2096    end do
2097    call xmpi_sum(resid_new(1),comm_atom_,ierr)
2098  end if
2099  vresid = resid_new(1)
2100 
2101 !_______________________________________________________________
2102 !Here use only the preconditioning, or initialize the other algorithms
2103 
2104  if(istep==1 .or. iscf==2)then
2105 
2106    write(message,'(2a)') ch10,'Simple mixing update:'
2107    call wrtout(std_out,message,'COLL')
2108 
2109    write(message,*)' residual square of the potential :',resid_new(1)
2110    call wrtout(std_out,message,'COLL')
2111 
2112 !  Store information for later use
2113    if (iscf==3.or.iscf==4) resid_old=resid_new(1)
2114    if (iscf==7) then
2115      amat(:,:)=zero
2116      amat(1,1)=resid_new(1)
2117    end if
2118 
2119 !  Compute new vtrial (and new rhoij if PAW)
2120    if (iscf/=2) f_fftgr(:,:,i_vstore)=vtrial(:,:)
2121    vtrial(:,:)=vtrial(:,:)+f_fftgr(:,:,i_vrespc(1))
2122    if (usepaw==1) then
2123      if (iscf/=2) f_paw(:,i_vstore)=vpaw(:)
2124      vpaw(:)=vpaw(:)+f_paw(:,i_vrespc(1))
2125    end if
2126 
2127 !  _______________________________________________________________
2128 !  Here Anderson algorithm using one previous iteration
2129  else if((istep==2 .or. iscf==3).and.iscf/=7)then
2130 
2131    write(message,'(2a)') ch10,'Anderson update:'
2132    call wrtout(std_out,message,'COLL')
2133 
2134    write(message,*)' residual square of the potential: ',resid_new(1)
2135    call wrtout(std_out,message,'COLL')
2136 
2137 !  Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2))
2138    call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,&
2139 &   nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2140    if (usepaw==1.and.pawoptmix==1) then
2141      do index=1,npawmix
2142        prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2))
2143      end do
2144      call xmpi_sum(prod_resid(1),comm_atom_,ierr)
2145    end if
2146 
2147 !  Compute mixing factor
2148    lambda=(resid_new(1)-prod_resid(1))/(resid_new(1)+resid_old-2*prod_resid(1))
2149    write(message,*)' mixing of old trial potential    :',lambda
2150    call wrtout(std_out,message,'COLL')
2151 
2152 !  Evaluate best residual square on the line
2153    resid_best=(1.0_dp-lambda)*(1.0_dp-lambda)*resid_new(1)&
2154 &   +(1.0_dp-lambda)*lambda        *2*prod_resid(1)&
2155 &   +lambda        *lambda        *resid_old
2156    write(message,*)' predicted best residual square on the line: ',resid_best
2157    call wrtout(std_out,message,'COLL')
2158 
2159 !  Store information for later use
2160    if (iscf==4) then
2161      prod_resid_old=prod_resid(1)
2162      resid_old2=resid_old
2163    end if
2164    resid_old=resid_new(1)
2165 
2166 !  Save latest trial potential and compute new trial potential
2167    do isp=1,nspden
2168      do ifft=1,cplex*nfft
2169        current=vtrial(ifft,isp)
2170        vtrial(ifft,isp)=(one-lambda)*(current                      +f_fftgr(ifft,isp,i_vrespc(1)))&
2171 &       +lambda      *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))
2172        f_fftgr(ifft,isp,i_vstore)=current
2173      end do
2174    end do
2175 
2176 !  PAW: save latest rhoij and compute new rhoij
2177    do index=1,npawmix
2178      current=vpaw(index)
2179      vpaw(index)=(one-lambda)*(current                 +f_paw(index,i_vrespc(1)))&
2180 &     +lambda      *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))
2181      f_paw(index,i_vstore)=current
2182    end do
2183 
2184 !  _______________________________________________________________
2185 !  Here Anderson algorithm using two previous iterations
2186  else if(iscf==4.and.iscf/=7)then
2187 
2188    write(message,'(2a)') ch10,'Anderson (order 2) update:'
2189    call wrtout(std_out,message,'COLL')
2190 
2191    write(message,*)' residual square of the potential :',resid_new(1)
2192    call wrtout(std_out,message,'COLL')
2193 
2194 !  Compute prod_resid from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(2))
2195    call dotprodm_v(cplex,1,prod_resid,i_vrespc(1),i_vrespc(2),mpicomm,mpi_summarize,1,1,&
2196 &   nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2197    if (usepaw==1.and.pawoptmix==1) then
2198      do index=1,npawmix
2199        prod_resid(1)=prod_resid(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(2))
2200      end do
2201    end if
2202 
2203 !  Compute prod_resid2 from f_fftgr/f_paw(:,:,i_vrespc(1)) and f_fftgr/f_paw(:,:,i_vrespc(3))
2204    call dotprodm_v(cplex,1,prod_resid2,i_vrespc(1),i_vrespc(3),mpicomm,mpi_summarize,1,1,&
2205 &   nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2206    if (usepaw==1.and.pawoptmix==1) then
2207      do index=1,npawmix
2208        prod_resid2(1)=prod_resid2(1)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(3))
2209      end do
2210 !    MPI reduction
2211      mpibuff(1)=prod_resid(1);mpibuff(2)=prod_resid2(1)
2212      call xmpi_sum(mpibuff,comm_atom_,ierr)
2213      prod_resid(1)=mpibuff(1);prod_resid2(1)=mpibuff(2)
2214    end if
2215 
2216 !  Compute mixing factors
2217    aa1=resid_new(1)+resid_old -two*prod_resid (1)
2218    aa2=resid_new(1)+resid_old2-two*prod_resid2(1)
2219    bb =resid_new(1)+prod_resid_old-prod_resid(1)-prod_resid2(1)
2220    cc1=resid_new(1)-prod_resid (1)
2221    cc2=resid_new(1)-prod_resid2(1)
2222    det=aa1*aa2-bb*bb
2223    lambda =(aa2*cc1-bb*cc2)/det
2224    lambda2=(aa1*cc2-bb*cc1)/det
2225    write(message,*)' mixing of old trial potentials   :',lambda,lambda2
2226    call wrtout(std_out,message,'COLL')
2227 
2228 !  Store information for later use
2229    prod_resid_old=prod_resid(1)
2230    resid_old2=resid_old
2231    resid_old=resid_new(1)
2232 
2233 !  Save latest trial potential and compute new trial potential
2234    do isp=1,nspden
2235      do ifft=1,cplex*nfft
2236        current=vtrial(ifft,isp)
2237        vtrial(ifft,isp)=&
2238 &       (one-lambda-lambda2)*(current                      +f_fftgr(ifft,isp,i_vrespc(1)))&
2239 &       +lambda             *(f_fftgr(ifft,isp,i_vtrial(1))+f_fftgr(ifft,isp,i_vrespc(2)))&
2240 &       +lambda2            *(f_fftgr(ifft,isp,i_vtrial(2))+f_fftgr(ifft,isp,i_vrespc(3)))
2241        f_fftgr(ifft,isp,i_vstore)=current
2242      end do
2243    end do
2244 
2245 !  PAW: save latest rhoij and compute new rhoij
2246    do index=1,npawmix
2247      current=vpaw(index)
2248      vpaw(index)=&
2249 &     (one-lambda-lambda2)*(current                 +f_paw(index,i_vrespc(1)))&
2250 &     +lambda             *(f_paw(index,i_vtrial(1))+f_paw(index,i_vrespc(2)))&
2251 &     +lambda2            *(f_paw(index,i_vtrial(2))+f_paw(index,i_vrespc(3)))
2252      f_paw(index,i_vstore)=current
2253    end do
2254 
2255 !  _______________________________________________________________
2256 !  Here Pulay algorithm
2257  else if(iscf==7)then
2258 
2259    niter=min(istep,npulay+1)
2260 
2261    write(message,'(2a,i2,a)') ch10,' Pulay update with ',niter-1,' previous iterations:'
2262    call wrtout(std_out,message,'COLL')
2263 
2264    if (npulay>npulaymax) then
2265      errid = AB7_ERROR_MIXING_CONVERGENCE
2266      write(errmess, '(4a)' ) ch10,&
2267 &     ' scfopt : ERROR - ',ch10,&
2268 &     '  Too much iterations required for Pulay algorithm (<50) !'
2269      return
2270    end if
2271 
2272 !  Compute "A" matrix
2273    if (istep>npulay+1) then
2274      do jj=1,niter-1
2275        do ii=1,niter-1
2276          amat(ii,jj)=amat(ii+1,jj+1)
2277        end do
2278      end do
2279    end if
2280    if (usepaw==1.and.pawoptmix==1) then
2281      ABI_ALLOCATE(amat_paw,(niter))
2282      amat_paw(:)=zero
2283      do ii=1,niter
2284        do index=1,npawmix
2285          amat_paw(ii)=amat_paw(ii)+f_paw(index,i_vrespc(1))*f_paw(index,i_vrespc(1+niter-ii))
2286        end do
2287      end do
2288      call xmpi_sum(amat_paw,comm_atom_,ierr)
2289    end if
2290    do ii=1,niter
2291      call dotprodm_v(cplex,1,amat(ii,niter),i_vrespc(1),i_vrespc(1+niter-ii),mpicomm,mpi_summarize,1,1,&
2292 &     nfft,n_fftgr,n_fftgr,nspden,opt_denpot,f_fftgr,f_fftgr)
2293      if (usepaw==1.and.pawoptmix==1) amat(ii,niter)=amat(ii,niter)+amat_paw(ii)
2294      if (ii<niter) amat(niter,ii)=amat(ii,niter)
2295    end do
2296    if (usepaw==1.and.pawoptmix==1)then
2297      ABI_DEALLOCATE(amat_paw)
2298    end if
2299 
2300 !  Invert "A" matrix
2301    ABI_ALLOCATE(amatinv,(niter,niter))
2302    amatinv(1:niter,1:niter)=amat(1:niter,1:niter)
2303    ABI_ALLOCATE(ipiv,(niter))
2304    ABI_ALLOCATE(rwork,(niter))
2305    call dgetrf(niter,niter,amatinv,niter,ipiv,ierr)
2306    call dgetri(niter,amatinv,niter,ipiv,rwork,niter,ierr)
2307    ABI_DEALLOCATE(ipiv)
2308    ABI_DEALLOCATE(rwork)
2309 
2310 !  Compute "alpha" factors
2311    ABI_ALLOCATE(alpha,(niter))
2312    alpha=zero
2313    det=zero
2314    do ii=1,niter
2315      do jj=1,niter
2316        alpha(ii)=alpha(ii)+amatinv(jj,ii)
2317        det=det+amatinv(jj,ii)
2318      end do
2319    end do
2320    alpha(:)=alpha(:)/det
2321    ABI_DEALLOCATE(amatinv)
2322    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)
2323    call wrtout(std_out,message,'COLL')
2324 
2325 !  Save latest trial potential and compute new trial potential
2326    do isp=1,nspden
2327      do ifft=1,cplex*nfft
2328        current=vtrial(ifft,isp)
2329        vtrial(ifft,isp)=alpha(niter)*(current+f_fftgr(ifft,isp,i_vrespc(1)))
2330        do ii=niter-1,1,-1
2331          vtrial(ifft,isp)=vtrial(ifft,isp)+alpha(ii) &
2332 &         *(f_fftgr(ifft,isp,i_vtrial(niter-ii))+f_fftgr(ifft,isp,i_vrespc(1+niter-ii)))
2333        end do
2334        f_fftgr(ifft,isp,i_vstore)=current
2335      end do
2336    end do
2337 
2338 !  PAW: save latest rhoij and compute new rhoij
2339    do index=1,npawmix
2340      current=vpaw(index)
2341      vpaw(index)=alpha(niter)*(current+f_paw(index,i_vrespc(1)))
2342      do ii=niter-1,1,-1
2343        vpaw(index)=vpaw(index)+alpha(ii) &
2344 &       *(f_paw(index,i_vtrial(niter-ii))+f_paw(index,i_vrespc(1+niter-ii)))
2345      end do
2346      f_paw(index,i_vstore)=current
2347    end do
2348 
2349    ABI_DEALLOCATE(alpha)
2350 !  _______________________________________________________________
2351 !  End of choice of optimization method
2352  end if
2353 
2354 !Permute potential indices
2355  if (iscf==3) then
2356    tmp=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
2357  else if (iscf==4) then
2358    tmp=i_vrespc(3) ; i_vrespc(3)=i_vrespc(2) ; i_vrespc(2)=i_vrespc(1) ; i_vrespc(1)=tmp
2359    tmp=i_vtrial(2) ; i_vtrial(2)=i_vtrial(1) ; i_vtrial(1)=tmp
2360  else if (iscf==7) then
2361    tmp=i_vtrial(  npulay)
2362    do ii=  npulay,2,-1
2363      i_vtrial(ii)=i_vtrial(ii-1)
2364    end do
2365    i_vtrial(1)=tmp
2366    tmp=i_vrespc(1+npulay)
2367    do ii=1+npulay,2,-1
2368      i_vrespc(ii)=i_vrespc(ii-1)
2369    end do
2370    i_vrespc(1)=tmp
2371  end if
2372 
2373 end subroutine scfopt