./72_response/dfpt_init

ABINIT/dfpt_init_mag1

ABINIT/dfpt_init_mag1 [ Functions ]

NAME

  dfpt_init_mag1

FUNCTION

  Initial guess of the first order magnetization/density for magnetic field perturbation.
  The first order magnetization is set so as to zero out the first order XC magnetic field, which
  should minimize the second order XC energy (without taking self-consistency into account).
  Works only for ipert==natom+5.

  Copyright (C) 2017-2018 ABINIT group (SPr)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

  idir  = direction of the applied magnetic field
  cplex = complex or real first order density and magnetization
  nfft  = dimension of the fft grid
  nspden= number of density matrix components
  nkxc  = number of kxc components
  vxc0(nfft,nspden)  = GS XC potential 
  kxc0(nfft,nspden)  = GS XC derivatives 
  rhor0(nfft,nspden) = GS density matrix

OUTPUT

  rhor1(cplex*nfft) = first order density magnetization guess

SIDE EFFECTS

NOTES

PARENTS

      dfpt_looppert

CHILDREN

SOURCE

 42 #if defined HAVE_CONFIG_H
 43 #include "config.h"
 44 #endif
 45 
 46 #include "abi_common.h"
 47 
 48 
 49 subroutine dfpt_init_mag1(idir,rhor1,rhor0,cplex,nfft,nspden,vxc0,kxc0,nkxc)
 50     
 51  use defs_basis
 52  use m_errors
 53  use m_profiling_abi
 54 
 55 !This section has been created automatically by the script Abilint (TD).
 56 !Do not modify the following lines by hand.
 57 #undef ABI_FUNC
 58 #define ABI_FUNC 'dfpt_init_mag1'
 59 !End of the abilint section
 60 
 61  implicit none
 62 
 63 !Arguments ------------------------------------
 64  integer , intent(in)    :: idir,cplex,nfft,nspden,nkxc
 65  real(dp), intent(in)    :: vxc0(nfft,nspden),rhor0(nfft,nspden)
 66  real(dp), intent(in)    :: kxc0(nfft,nkxc)
 67  real(dp), intent(out)   :: rhor1(cplex*nfft,nspden)                        
 68 
 69 !Local variables-------------------------------
 70  integer  :: ipt                                     
 71  real(dp) :: bxc0,bxc1                  
 72  real(dp) :: m1_norm,m0_norm
 73  real(dp) :: f_dot_m
 74  real(dp) :: mdir(3),fdir(3)               
 75  
 76 ! *************************************************************************
 77 
 78  if (nspden==2) then
 79 
 80    if(cplex==1) then
 81      do ipt=1,nfft
 82        bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2)) ! d/dm Bxc
 83        !this overestimates the first order magnetization because of n1 not taken into account
 84        m1_norm=-half*(1/bxc1)
 85        rhor1(ipt,1)=zero             ! rho_up+rho_dwn    => charge density
 86        rhor1(ipt,2)=half*m1_norm     ! rho_up=1/2(rho+m) => half*m
 87      end do
 88    else
 89      do ipt=1,cplex*nfft
 90        rhor1(ipt,:)=zero
 91      end do
 92    end if
 93 
 94  else if(nspden==4) then
 95 
 96    fdir=zero
 97    fdir(idir)= 1.0d0
 98    do ipt=1,nfft  
 99      m0_norm=sqrt(rhor0(ipt,2)**2+rhor0(ipt,3)**2+rhor0(ipt,4)**2)
100      mdir(1)=rhor0(ipt,2)/m0_norm
101      mdir(2)=rhor0(ipt,3)/m0_norm
102      mdir(3)=rhor0(ipt,4)/m0_norm
103      f_dot_m=fdir(1)*mdir(1)+fdir(2)*mdir(2)+fdir(3)*mdir(3) ! projection of the field direction on m0
104 
105      bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2))  ! d/dm Bxc
106      m1_norm=(-half/bxc1)*f_dot_m                            ! get an estimate of the norm of m1
107 
108      bxc0=-sqrt((half*(vxc0(ipt,1)-vxc0(ipt,2)))**2+vxc0(ipt,3)**2+vxc0(ipt,4)**2)       
109 
110      if(cplex==1) then
111        rhor1(ipt,1)=zero       ! rho_up+rho_dwn    => charge density
112        rhor1(ipt,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1))   ! m1x
113        rhor1(ipt,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2))   ! m1x
114        rhor1(ipt,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3))   ! m1x
115      else
116        rhor1(2*ipt-1,1)=zero       ! Re rho_up+rho_dwn
117        rhor1(2*ipt-1,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1))   ! m1x
118        rhor1(2*ipt-1,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2))   ! m1x
119        rhor1(2*ipt-1,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3))   ! m1x
120        rhor1(2*ipt  ,1)=zero       ! Im rho_up+rho_dwn
121        rhor1(2*ipt  ,2)=zero
122        rhor1(2*ipt  ,3)=zero
123        rhor1(2*ipt  ,4)=zero
124 
125        rhor1(2*ipt-1,1)=zero; rhor1(2*ipt,1)=zero
126        rhor1(2*ipt-1,2)=zero; rhor1(2*ipt,2)=zero
127        rhor1(2*ipt-1,3)=zero; rhor1(2*ipt,3)=zero
128        rhor1(2*ipt-1,4)=zero; rhor1(2*ipt,4)=zero
129 
130      end if
131    end do
132  end if
133 
134 
135 end subroutine dfpt_init_mag1

TABLE OF CONTENTS

ABINIT/dfpt_init_mag1 [ Functions ]