TABLE OF CONTENTS


ABINIT/m_rhotov [ Modules ]

[ Top ] [ Modules ]

NAME

  m_rhotov

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (XG, GMR, MT, EB)
  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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_rhotov
27 
28  use defs_basis
29  use defs_abitypes
30  use defs_wvltypes
31  use m_errors
32  use m_abicore
33  use m_ab7_mixing
34  use m_abi2big
35  use m_xmpi
36  use m_cgtools
37  use m_xcdata
38 
39  use m_time,             only : timab
40  use m_geometry,         only : xred2xcart
41  use m_energies,         only : energies_type
42  use m_electronpositron, only : electronpositron_type, electronpositron_calctype, rhohxcpositron
43  use libxc_functionals,  only : libxc_functionals_is_hybrid
44  use m_spacepar,         only : hartre
45  use m_dens,             only : mag_constr
46  use m_rhotoxc,          only : rhotoxc
47  use m_xchybrid,         only : xchybrid_ncpp_cc
48  use m_psolver,          only : psolver_rhohxc
49  use m_wvl_psi,          only : wvl_psitohpsi
50 
51  implicit none
52 
53  private

ABINIT/rhotov [ Functions ]

[ Top ] [ Functions ]

NAME

 rhotov

FUNCTION

 This routine is called to compute, from a given total density
 the trial (local) potential and the residual potential.

INPUTS

  [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
  dtset <type(dataset_type)>=all input variables in this dataset
   | spinmagntarget=input variable that governs fixed moment calculation
   | natom=number of atoms in cell.
   | nspden=number of spin-density components
   | ntypat=number of types of atoms in unit cell.
   | occopt=option for occupancies
   | typat(natom)=type (integer) for each atom
  gprimd(3,3)=dimensional reciprocal space primitive translations
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  mpi_enreg=informations about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhat(nfft,nspden*usepaw)= -PAW only- compensation density
  nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
  nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
  nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description
  n3xccc=dimension of the xccc3d array (0 or nfft).
  optene=option for the computation of additional energies
  optres=0: the trial potential residual is computed ; the input potential value is kept
         1: the new value of the trial potential is computed in place of the input value
  optxc=option to be used for the call to rhotoxc
  rhog(2,nfft)=array for Fourier transform of electron density
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
   | definition for spin components:
   | case of nspden = 2
   |      rhor(:,1) => rho_up + rho_dwn
   |      rhor(:,2) => rho_up
   | case of nspden = 4
   |      rhor(:,1)   => rho_upup + rho_dwndwn
   |      rhor(:,2:4) => {m_x,m_y,m_z}
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  [taug(2,nfftf*dtset%usekden)]=array for Fourier transform of kinetic energy density
  [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density
  ucvol = unit cell volume (Bohr**3)
  usepaw= 0 for non paw calculation; =1 for paw calculation
  usexcnhat= -PAW only- flag controling use of compensation density in Vxc
  vpsp(nfft)=array for holding local psp
  [vxc_hybcomp(nfft,nspden)= compensation xc potential (Hartree) in case of hybrids] Optional output
       i.e. difference between the hybrid Vxc at fixed density and the auxiliary Vxc at fixed density
  xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
  ==== if optres==0
    vtrial(nfft,nspden)= old value of trial potential

OUTPUT

  energies <type(energies_type)>=all part of total energy.
   | e_hartree=Hartree part of total energy (hartree units)
   | e_xc=exchange-correlation energy (hartree)
   | In case of hybrid compensation algorithm:
   | e_hybcomp_v=self-consistent potential compensation term for the exchange-correlation energy (hartree)
  ==== if optene==0.or.2
   | e_localpsp=local psp energy (hartree)
  ==== if optene==1.or.2
   | e_xcdc=exchange-correlation double-counting energy (hartree)
  kxc(nfft,nkxc)=exchange-correlation kernel, needed only if optxc==2.
  strsxc(6)=xc contribution to stress tensor (hartree/bohr^3)
  vxc(nfft,nspden)=Vxc(r) (already computed above; gets recomputed below too)
  vxcavg=mean of the vxc potential
  ==== if optres==0
    vresidnew(nfft,nspden)=potential residual
    vnew_mean(nspden)=mean of the potential formed from vpsp, vhartr and vxc, might be spin-dependent
    vres_mean(nspden)=mean of the potential residual, might be spin-dependent
    vres2=square of the norm of the residual
    [vxctau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density with respect to
      kinetic energy density (metaGGA cases) (optional output)

SIDE EFFECTS

 Input/Output:
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
  vhartr(nfft)=array for holding Hartree potential
  ==== if optres==1
    vtrial(nfft,nspden)= new value of trial potential

NOTES

  In case of PAW calculations:
    All computations are done on the fine FFT grid.
    All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
    All arrays (densities/potentials...) are computed on this fine FFT grid.
  ! Developpers have to be careful when introducing others arrays:
      they have to be stored on the fine FFT grid.
  In case of norm-conserving calculations the FFT grid is the usual FFT grid.

PARENTS

      scfcv

CHILDREN

      dotprod_vn,hartre,mag_constr,mean_fftr,psolver_rhohxc,rhohxcpositron
      rhotoxc,sqnorm_v,timab,wvl_psitohpsi,wvl_vtrial_abi2big,xcdata_init
      xchybrid_ncpp_cc,xred2xcart

SOURCE

163 subroutine rhotov(dtset,energies,gprimd,gsqcut,istep,kxc,mpi_enreg,nfft,ngfft,&
164 &  nhat,nhatgr,nhatgrdim,nkxc,vresidnew,n3xccc,optene,optres,optxc,&
165 &  rhog,rhor,rprimd,strsxc,ucvol,usepaw,usexcnhat,&
166 &  vhartr,vnew_mean,vpsp,vres_mean,vres2,vtrial,vxcavg,vxc,wvl,xccc3d,xred,&
167 &  electronpositron,taug,taur,vxc_hybcomp,vxctau,add_tfw) ! optional arguments
168 
169 
170 !This section has been created automatically by the script Abilint (TD).
171 !Do not modify the following lines by hand.
172 #undef ABI_FUNC
173 #define ABI_FUNC 'rhotov'
174 !End of the abilint section
175 
176  implicit none
177 
178 !Arguments ------------------------------------
179 !scalars
180  integer,intent(in) :: n3xccc,nfft,nhatgrdim,nkxc,optene,optres,optxc,usepaw,istep
181  integer,intent(in) :: usexcnhat
182  logical,intent(in),optional :: add_tfw
183  real(dp),intent(in) :: gsqcut,ucvol
184  real(dp),intent(out) :: vres2,vxcavg
185  type(MPI_type),intent(inout) :: mpi_enreg
186  type(dataset_type),intent(in) :: dtset
187  type(electronpositron_type),pointer,optional :: electronpositron
188  type(energies_type),intent(inout) :: energies
189  type(wvl_data), intent(inout) :: wvl
190 !arrays
191  integer,intent(in) :: ngfft(18)
192  real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*usepaw)
193  real(dp),intent(in) :: nhatgr(nfft,dtset%nspden,3*nhatgrdim),rhog(2,nfft)
194  real(dp),intent(in) :: rprimd(3,3)
195  real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft)
196  real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden)
197  real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom)
198  real(dp),intent(out) :: kxc(nfft,nkxc),strsxc(6),vnew_mean(dtset%nspden)
199  real(dp),intent(out) :: vres_mean(dtset%nspden),vresidnew(nfft,dtset%nspden)
200  real(dp),intent(in),optional :: taug(2,nfft*dtset%usekden)
201  real(dp),intent(in),optional :: taur(nfft,dtset%nspden*dtset%usekden)
202  real(dp),intent(out),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden)
203  real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden)
204 
205 !Local variables-------------------------------
206 !scalars
207  integer :: nk3xc,ifft,ipositron,ispden,nfftot,offset
208  integer :: mpi_comm_sphgrid,ixc_current
209 !integer :: ii,jj,kk,ipt,nx,ny,nz           !SPr: debug
210 !real(dp):: rx,ry,rz                        !SPr: debug
211  real(dp) :: doti,e_xcdc_vxctau
212  logical :: add_tfw_,calc_xcdc,non_magnetic_xc,with_vxctau
213  logical :: is_hybrid_ncpp,wvlbigdft=.false.
214  type(xcdata_type) :: xcdata
215 !arrays
216  real(dp) :: evxc,tsec(2),vmean(dtset%nspden),vzeeman(dtset%nspden)
217  real(dp),allocatable :: rhowk(:,:),Vmagconstr(:,:),vnew(:,:),xcart(:,:)
218 !real(dp),allocatable :: vzeemanHarm(:,:)   !SPr: debug Zeeman field q/=0 real space
219 
220 ! *********************************************************************
221 
222  DBG_ENTER("COLL")
223 
224  call timab(940,1,tsec)
225 
226 ! Initialise non_magnetic_xc for rhohxc
227  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
228 
229 !Check that usekden is not 0 if want to use vxctau
230  with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0))
231 
232 !Check if we're in hybrid norm conserving pseudopotential with a core correction
233  is_hybrid_ncpp=(usepaw==0 .and. n3xccc/=0 .and. &
234 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
235 
236 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
237  wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
238 
239 !mpi communicator for spherical grid
240  mpi_comm_sphgrid=mpi_enreg%comm_fft
241  if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
242 
243 !Get size of FFT grid
244  nfftot=PRODUCT(ngfft(1:3))
245 
246  ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron)
247  add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
248 
249 !------Compute Hartree and xc potentials----------------------------------
250 
251 !allocate vnew here.
252 !In wvl: vnew is used at call to wvl_psitohpsi
253  if (optres==0) then
254    ABI_ALLOCATE(vnew,(nfft,dtset%nspden))
255    vmean(:)=zero ; vnew_mean(:)=zero
256  end if
257 
258  if (ipositron/=1) then
259 !  Compute xc potential (separate up and down if spin-polarized)
260    if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then
261      call hartre(1,gsqcut,usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartr)
262      !Use the proper exchange_correlation energy : either the origin one, or the auxiliary one
263      ixc_current=dtset%ixc
264      if(mod(dtset%fockoptmix,100)==11)ixc_current=dtset%auxc_ixc
265      call xcdata_init(xcdata,dtset=dtset,ixc=ixc_current)
266 
267 !    Use the periodic solver to compute Hxc.
268      nk3xc=1
269 !write(80,*) "rhotov"
270 !xccc3d=zero
271 !DEBUG
272      write(std_out,*)' rhotov : is_hybrid_ncpp=',is_hybrid_ncpp
273      write(std_out,*)' rhotov : n3xccc=',n3xccc
274 !ENDDEBUG
275 
276      call timab(941,1,tsec)
277      if (ipositron==0) then
278        if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then
279          call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,&
280 &         nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,dtset%paral_kgb,&
281 &         rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,&
282 &         taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_)
283          if(mod(dtset%fockoptmix,100)==11)then
284            energies%e_xc=energies%e_xc*dtset%auxc_scal
285            vxc(:,:)=vxc(:,:)*dtset%auxc_scal
286          end if
287        else
288          call xchybrid_ncpp_cc(dtset,energies%e_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
289 &         strsxc,vxcavg,xccc3d,vxc=vxc)
290        end if
291      else
292        call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,&
293 &       nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,dtset%paral_kgb,&
294 &       rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,&
295 &       taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,&
296 &       electronpositron=electronpositron)
297      end if
298      call timab(941,2,tsec)
299 
300    elseif (.not. wvlbigdft) then
301 !    Use the free boundary solver.
302      call timab(943,1,tsec)
303      call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
304 &     dtset%icoulomb, dtset%ixc, &
305 &     mpi_enreg, nfft, &
306 &     ngfft, nhat,usepaw,&
307 &     dtset%nscforder, dtset%nspden, n3xccc, rhor,rprimd,&
308 &     usexcnhat,dtset%usepaw,dtset%usewvl,vhartr, vxc, vxcavg,&
309 &     wvl%descr,wvl%den,wvl%e,&
310 &     xccc3d,dtset%xclevel,dtset%xc_denpos)
311      call timab(943,2,tsec)
312    end if
313 !  For icoulomb==0 and usewvl Ehartree is calculated in psolver_rhohxc().
314 !  For PAW we recalculate this since nhat was not taken into account
315 !  in psolver_rhohxc: E_H= int v_H (n+nhat) dr
316    if(.not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then
317 
318      call timab(942,1,tsec)
319      call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
320      energies%e_hartree=half*energies%e_hartree
321      call timab(942,2,tsec)
322    end if
323  else
324    call timab(944,1,tsec)
325    energies%e_hartree=zero;energies%e_xc=zero
326    call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,&
327 &   dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos)
328    call timab(944,2,tsec)
329  end if
330 
331  call timab(945,1,tsec)
332  if (ipositron/=0) then
333    call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,&
334 &   nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
335    vhartr=vhartr+electronpositron%vha_ep
336  end if
337 
338 !------Compute parts of total energy depending on potentials--------
339 
340  if ( (optene==0.or.optene==2 ).and. .not. wvlbigdft) then
341 !  Compute local psp energy energies%e_localpsp
342    call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol,&
343 &   mpi_comm_sphgrid=mpi_comm_sphgrid)
344  end if
345 
346  if(mod(dtset%fockoptmix,100)==11)then
347    if (.not. wvlbigdft) then
348 !    Compute second compensation energy for hybrid functionals
349      call dotprod_vn(1,rhor,energies%e_hybcomp_v,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol,&
350 &     mpi_comm_sphgrid=mpi_comm_sphgrid)
351    end if
352  end if
353 
354  calc_xcdc=.false.
355  if (optene==1.or.optene==2) calc_xcdc=.true.
356  if (dtset%usewvl==1.and.dtset%nnsclo>0) calc_xcdc=.true.
357  if (wvlbigdft) calc_xcdc=.false.
358  if (dtset%usefock==1) calc_xcdc=.true.
359 
360  if (calc_xcdc) then
361 
362 !  Compute double-counting XC energy energies%e_xcdc
363    if (ipositron/=1) then
364      if (usepaw==0.or.usexcnhat/=0) then
365        call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,&
366 &       mpi_comm_sphgrid=mpi_comm_sphgrid)
367        if (with_vxctau)then
368          call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),&
369 &         ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
370          energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
371        end if
372      else
373        ABI_ALLOCATE(rhowk,(nfft,dtset%nspden))
374        rhowk=rhor-nhat
375        call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,&
376 &       mpi_comm_sphgrid=mpi_comm_sphgrid)
377        ABI_DEALLOCATE(rhowk)
378      end if
379      if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc
380    else
381      energies%e_xcdc=zero
382    end if
383 
384  end if
385 
386 !------Produce residual vector and square norm of it-------------
387 !(only if requested ; if optres==0)
388 
389 !Set up array for Zeeman field
390 !EB vzeeman(:) = factor*( Hz, Hx+iHy; Hx-iHy, -Hz)
391 !EB factor = -g/2 * mu_B * mu_0 = -1/2*B in a.u.
392 !EB <-- vzeeman might have to be allocated correctly --> to be checked
393 ! vzeeman = 1/2 ( -B_z, -B_x + iB_y ; -B_x - iB_y , B_z)
394  vzeeman(:) = zero
395 ! ABI_ALLOCATE(vzeemanHarm,(nfft,dtset%nspden))  ! SPr: debug stuff
396 ! vzeemanHarm(:,:) = zero                        !
397  if (any(abs(dtset%zeemanfield(:))>tol8)) then
398    if(dtset%nspden==2)then
399 !    EB The collinear case has to be checked :
400 !    EB Is it vzeeman(1) or (2) that has to be added here? to be checked in setvtr and energy as well
401 !    SPr: the density components are: rhor(1) => n_upup + n_dwndwn
402 !                                     rhor(2) => n_upup
403 !         the convention for the potential components is different:
404 !                                     v(1)    => v_upup
405 !                                     v(2)    => v_dndn
406 !         verified by comparing collinear and non-collinear calculations
407 
408      vzeeman(1) =-half*dtset%zeemanfield(3)  ! v_upup
409      vzeeman(2) = half*dtset%zeemanfield(3)  ! v_dndn
410 
411      !vzeeman(1) = zero  ! v_upup
412      !vzeeman(2) = zero  ! v_dndn
413 
414      !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3)
415      !do kk=0,nz-1
416      !  do jj=0,ny-1
417      !    do ii=0,nx-1
418      !      ipt=1+ii+nx*(jj+ny*kk)
419      !      !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3)
420      !      !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3)
421      !      !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3)
422      !      vzeemanHarm(ipt,1)= -half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
423      !      vzeemanHarm(ipt,2)=  half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
424      !    end do
425      !  end do
426      !end do
427 
428    else if(dtset%nspden==4)then
429 
430      vzeeman(1)=-half*dtset%zeemanfield(3)    ! v_upup
431      vzeeman(2)= half*dtset%zeemanfield(3)    ! v_dndn
432      vzeeman(3)=-half*dtset%zeemanfield(1)    ! Re(v_updn)
433      vzeeman(4)= half*dtset%zeemanfield(2)    ! Im(v_updn)
434 
435      !vzeeman(1)=0.0
436      !vzeeman(2)=0.0
437      !vzeeman(3)=0.0
438      !vzeeman(4)=0.0
439 
440      !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3)
441      !do kk=0,nz-1
442      !  do jj=0,ny-1
443      !    do ii=0,nx-1
444      !      ipt=1+ii+nx*(jj+ny*kk)
445      !      !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3)
446      !      !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3)
447      !      !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3)
448      !      vzeemanHarm(ipt,1)= -half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
449      !      vzeemanHarm(ipt,2)=  half*dtset%zeemanfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
450      !      vzeemanHarm(ipt,3)= -half*dtset%zeemanfield(1)*cos(2*PI*(dble(ii)/dble(nx)))
451      !      vzeemanHarm(ipt,4)=  half*dtset%zeemanfield(2)*cos(2*PI*(dble(ii)/dble(nx)))
452      !    end do
453      !  end do
454      !end do
455 
456    end if
457  end if
458 
459 !Compute the constrained potential for the magnetic moments
460  ABI_ALLOCATE(Vmagconstr, (nfft,dtset%nspden))
461  if (dtset%magconon==1.or.dtset%magconon==2) then
462    Vmagconstr = zero
463    call mag_constr(dtset%natom,dtset%spinat,dtset%nspden,dtset%magconon,dtset%magcon_lambda,rprimd, &
464 &   mpi_enreg,nfft,ngfft,dtset%ntypat,dtset%ratsph,rhor,dtset%typat,Vmagconstr,xred)
465  else
466 !  NOTE: mjv 25 May 2013 added this for ibm6 - otherwise gives NaN in vnew after
467 !  the addition below, in case magconon==0 and for certain libxc
468 !  functionals!! May need to copy this to setvtr.F90 if the same effect appears.
469 !  should check libxc test with a proper memory checker (valgrind).
470    do ispden=1,dtset%nspden
471      do ifft=1,nfft
472        Vmagconstr(ifft,ispden) = zero
473      end do
474    end do
475  end if
476 
477  if (optres==0) then
478 
479 
480 !  ------ Compute potential residual -------------
481 
482    if (.not. wvlbigdft) then
483 !$OMP PARALLEL DO COLLAPSE(2)
484      do ispden=1,min(dtset%nspden,2)
485        do ifft=1,nfft
486          vnew(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vzeeman(ispden)+Vmagconstr(ifft,ispden)
487          !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden)
488          if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden)
489          vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden)
490        end do
491      end do
492      if(dtset%nspden==4)then
493 !$OMP PARALLEL DO COLLAPSE(2)
494        do ispden=3,4
495          do ifft=1,nfft
496            vnew(ifft,ispden)=vxc(ifft,ispden)+vzeeman(ispden)+Vmagconstr(ifft,ispden)
497            !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden)
498            if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden)
499            vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden)
500          end do
501        end do
502      end if
503 
504      offset   = 0
505 
506      if (dtset%iscf==0) vtrial=vnew
507 
508 !    Pass vtrial to BigDFT object
509      if(dtset%usewvl==1) then
510        call wvl_vtrial_abi2big(1,vnew,wvl%den)
511 !      call wvl_vtrial_abi2big(1,vtrial,wvl%den)
512      end if
513 
514    else
515 !    Compute with covering comms the different part of the potential.
516 !    only for wvlbigdft
517      ABI_ALLOCATE(xcart,(3, dtset%natom))
518      call xred2xcart(dtset%natom, rprimd, xcart, xred)
519      call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, &
520 &     energies%e_kinetic, energies%e_localpsp, energies%e_nonlocalpsp, energies%e_sicdc, &
521 &     istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft,&
522 &     mpi_enreg%nproc_wvl, dtset%nspden, &
523 &     vres2, .true., energies%e_xcdc, wvl,&
524 &     wvlbigdft, xcart, strsxc,vtrial=vnew,vxc=vxc)
525      ABI_DEALLOCATE(xcart)
526 
527      vresidnew = vnew - vtrial
528      vtrial = vnew
529 
530      call mean_fftr(vxc, vmean(1:1),  nfft, nfftot, dtset%nspden,&
531 &     mpi_comm_sphgrid=mpi_comm_sphgrid)
532      vxcavg = vmean(1)
533      offset = 0
534    end if
535 
536 !  Compute mean values of potential and residual
537    call mean_fftr(vnew(1+offset, 1),vnew_mean,nfft,nfftot,dtset%nspden,&
538 &   mpi_comm_sphgrid=mpi_comm_sphgrid)
539    call mean_fftr(vresidnew(1+offset, 1),vmean,nfft,nfftot,dtset%nspden,&
540 &   mpi_comm_sphgrid=mpi_comm_sphgrid)
541 
542    ABI_DEALLOCATE(vnew)
543 
544 !  Subtract the mean of the residual
545 !  Must take into account fixed occupation number in case of spin-polarized
546    do ispden=1,dtset%nspden
547      if (dtset%nspden==2.and.dtset%occopt>=3.and. abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then
548        vres_mean(ispden)=(vmean(1)+vmean(2))*half
549      else
550        vres_mean(ispden)=vmean(ispden)
551      end if
552    end do
553 
554 !$OMP PARALLEL DO COLLAPSE(2)
555    do ispden=1,dtset%nspden
556      do ifft=1,nfft
557        vresidnew(ifft,ispden)=vresidnew(ifft,ispden)-vres_mean(ispden)
558      end do
559    end do
560 
561 !  Compute square norm vres2 of potential residual vresid
562    call sqnorm_v(1,nfft,vres2,dtset%nspden,optres,vresidnew(1+offset, 1),mpi_comm_sphgrid=mpi_comm_sphgrid)
563 
564  else ! optres/=0
565 
566 !  ------Produce new value of trial potential-------------
567 
568    if (.not. wvlbigdft) then
569 !$OMP PARALLEL DO COLLAPSE(2)
570      do ispden=1,min(dtset%nspden,2)
571        do ifft=1,nfft
572          vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vzeeman(ispden)+Vmagconstr(ifft,ispden)
573          !vtrial(ifft,ispden)=vtrial(ifft,ispden)+vzeemanHarm(ifft,ispden)
574          if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,ispden)=vtrial(ifft,ispden)+vxc_hybcomp(ifft,ispden)
575        end do
576      end do
577      if(dtset%nspden==4) then
578 !$OMP PARALLEL DO
579        do ifft=1,nfft
580          vtrial(ifft,3:4)=vxc(ifft,3:4)+vzeeman(3:4)+Vmagconstr(ifft,3:4)
581          !vtrial(ifft,3:4)=vtrial(ifft,3:4)+vzeemanHarm(ifft,3:4)
582          if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,3:4)=vtrial(ifft,3:4)+vxc_hybcomp(ifft,3:4)
583        end do
584      end if
585 !    Pass vtrial to BigDFT object
586      if(dtset%usewvl==1) then
587        call wvl_vtrial_abi2big(1,vtrial,wvl%den)
588      end if
589    else
590 !    Compute with covering comms the different part of the potential.
591      ABI_ALLOCATE(xcart,(3, dtset%natom))
592      call xred2xcart(dtset%natom, rprimd, xcart, xred)
593      call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, &
594 &     energies%e_kinetic, energies%e_localpsp, energies%e_nonlocalpsp, energies%e_sicdc, &
595 &     istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, &
596 &     dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl,&
597 &     dtset%nspden,vres2, .true.,energies%e_xcdc,  wvl,&
598 &     wvlbigdft, xcart, strsxc, vtrial, vxc)
599      ABI_DEALLOCATE(xcart)
600 !    Compute vxcavg
601      call mean_fftr(vxc, vmean(1:1), nfft, nfftot, dtset%nspden,&
602 &     mpi_comm_sphgrid=mpi_comm_sphgrid)
603      vxcavg = vmean(1)
604    end if
605 
606  end if
607 
608  ABI_DEALLOCATE(Vmagconstr)
609  !ABI_DEALLOCATE(vzeemanHarm) !SPr: debug for q/=0 magnetic field
610 
611  call timab(945,2,tsec)
612  call timab(940,2,tsec)
613 
614  DBG_EXIT("COLL")
615 
616 end subroutine rhotov