## 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
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
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
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
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',&
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',&
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