TABLE OF CONTENTS
ABINIT/mag_constr [ Functions ]
NAME
mag_constr
FUNCTION
This routine is called to compute the potential corresponding to constrained magnetic moments.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (ILuk, MVer) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
natom=number of atoms spinat=fixed magnetic moments vectors nspden = number of spin densities (1 2 or 4) magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size magcon_lambda=the size of the penalty terms rprimd=lattice vectors (dimensionful) mpi_enreg=mpi structure with communicator info nfft=number of points in standard fft grid ngfft=FFT grid dimensions ntypat=number of types of atoms ratsph=radii for muffin tin spheres of each atom rhor=density in real space typat=types of atoms xred=reduced atomic positions
OUTPUT
Vmagconstr=the constraining potential
PARENTS
energy,rhotov,setvtr
CHILDREN
calcdensph,metric,ptabs_fourdp,timab,xmpi_sum
NOTES
based on html notes for the VASP implementation at http://cms.mpi.univie.ac.at/vasp/vasp/Constraining_direction_magnetic_moments.html
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 subroutine mag_constr(natom,spinat,nspden,magconon,magcon_lambda,rprimd, & 54 mpi_enreg,nfft,ngfft,ntypat,ratsph,rhor, & 55 typat,Vmagconstr,xred) 56 57 use defs_basis 58 use defs_abitypes 59 use m_xmpi 60 use m_errors 61 62 use m_geometry, only : metric 63 use m_mpinfo, only : ptabs_fourdp 64 65 !This section has been created automatically by the script Abilint (TD). 66 !Do not modify the following lines by hand. 67 #undef ABI_FUNC 68 #define ABI_FUNC 'mag_constr' 69 use interfaces_18_timing 70 use interfaces_32_util 71 use interfaces_54_abiutil 72 !End of the abilint section 73 74 implicit none 75 76 !Arguments ------------------------------------ 77 !scalars 78 integer,intent(in) :: natom,magconon,nfft,nspden 79 integer,intent(in) :: ntypat 80 real(dp),intent(in) :: magcon_lambda 81 real(dp),intent(out) :: Vmagconstr(nfft,nspden) 82 type(MPI_type),intent(in) :: mpi_enreg 83 !arrays 84 integer,intent(in) :: typat(natom) 85 integer,intent(in) :: ngfft(18) 86 real(dp),intent(in) :: ratsph(ntypat) 87 real(dp),intent(in) :: rhor(nfft,nspden) 88 real(dp),intent(in) :: rprimd(3,3) 89 real(dp),intent(in) :: spinat(3,natom) 90 real(dp),intent(in) :: xred(3,natom) 91 92 !Local variables------------------------------- 93 !scalars 94 integer,parameter :: ishift=5 95 integer :: iatom, ierr 96 integer :: cplex1=1 97 integer :: n1a, n1b, n3a, n3b, n2a, n2b 98 integer :: n1, n2, n3 99 integer :: ifft_local 100 integer :: i1,i2,i3,ix,iy,iz,izloc,nd3 101 real(dp) :: arg,intgden_proj,r2atsph,rr1,rr2,rr3,fsm,ratsm,ratsm2,difx,dify,difz,r2,rx,ry,rz 102 real(dp) :: norm 103 real(dp),parameter :: delta=0.99_dp 104 !arrays 105 real(dp) :: intgden(nspden,natom) 106 real(dp) :: spinat_norm(3,natom) 107 !real(dp),allocatable :: cmm_x(:,:),cmm_y(:,:),cmm_z(:,:) 108 real(dp):: cmm_x,cmm_y,cmm_z 109 real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3),ucvol 110 real(dp) :: tsec(2) 111 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 112 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 113 114 ! *********************************************************************************************** 115 116 !We need the metric because it is needed in calcdensph.F90 117 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 118 119 !We need the integrated magnetic moments and the smoothing function 120 call calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,std_out,ratsph,rhor,rprimd,typat,ucvol,xred,1,cplex1,intgden=intgden) 121 122 n1 = ngfft(1) 123 n2 = ngfft(2) 124 n3 = ngfft(3) 125 nd3=n3/mpi_enreg%nproc_fft 126 127 ratsm = 0.05_dp ! default value for the smearing region radius - may become input variable later 128 Vmagconstr = zero 129 130 !Get the distrib associated with this fft_grid 131 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 132 133 !Loop over atoms 134 !------------------------------------------- 135 do iatom=1,natom 136 137 norm = sqrt(sum(spinat(:,iatom)**2)) 138 spinat_norm(:,iatom) = zero 139 if (norm > tol10) then 140 spinat_norm(:,iatom) = spinat(:,iatom) / norm 141 else if (magconon == 1) then 142 ! if spinat = 0 and we are imposing the direction only, skip this atom 143 cycle 144 end if 145 146 ! Define a "box" around the atom 147 r2atsph=1.0000001_dp*ratsph(typat(iatom))**2 148 rr1=sqrt(r2atsph*gmet(1,1)) 149 rr2=sqrt(r2atsph*gmet(2,2)) 150 rr3=sqrt(r2atsph*gmet(3,3)) 151 152 n1a=int((xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1 153 n1b=int((xred(1,iatom)+rr1+ishift)*n1 )-ishift*n1 154 n2a=int((xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2 155 n2b=int((xred(2,iatom)+rr2+ishift)*n2 )-ishift*n2 156 n3a=int((xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3 157 n3b=int((xred(3,iatom)+rr3+ishift)*n3 )-ishift*n3 158 159 ratsm2 = -(ratsm**2 - 2*ratsph(typat(iatom))*ratsm) 160 161 do i3=n3a,n3b 162 iz=mod(i3+ishift*n3,n3) 163 if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then 164 izloc = ffti3_local(iz+1) - 1 165 difz=dble(i3)/dble(n3)-xred(3,iatom) 166 do i2=n2a,n2b 167 iy=mod(i2+ishift*n2,n2) 168 dify=dble(i2)/dble(n2)-xred(2,iatom) 169 do i1=n1a,n1b 170 ix=mod(i1+ishift*n1,n1) 171 difx=dble(i1)/dble(n1)-xred(1,iatom) 172 rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3) 173 ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3) 174 rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3) 175 r2=rx**2+ry**2+rz**2 176 177 178 ! Identify the fft indexes of the rectangular grid around the atom 179 if (r2 > r2atsph) cycle 180 181 fsm = radsmear(r2, r2atsph, ratsm2) 182 183 ifft_local=1+ix+n1*(iy+n2*izloc) 184 185 ! Calculate the x- and y-components of the square bracket term 186 cmm_x = zero 187 cmm_y = zero 188 cmm_z = zero 189 intgden_proj = zero 190 if (nspden == 4) then 191 if (magconon==1) then 192 ! Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector 193 ! This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector 194 intgden_proj=spinat_norm(1,iatom)*intgden(2,iatom)+ & 195 & spinat_norm(2,iatom)*intgden(3,iatom)+ & 196 & spinat_norm(3,iatom)*intgden(4,iatom) 197 198 cmm_x=intgden(2,iatom) 199 cmm_x=cmm_x-spinat_norm(1,iatom)*intgden_proj 200 201 cmm_y=intgden(3,iatom) 202 cmm_y=cmm_y-spinat_norm(2,iatom)*intgden_proj 203 204 else if (magconon==2 .and. nspden == 4) then 205 cmm_x=intgden(2,iatom)-spinat(1,iatom) 206 cmm_y=intgden(3,iatom)-spinat(2,iatom) 207 end if 208 209 ! do loop on spirals - should be added later on. arg = scalar product k.R (with 2pi factor?) 210 arg=0 211 212 ! Calculate the constraining potential for x- and y- components of the mag. mom. vector 213 ! Eric Bousquet has derived the relationship between spin components and potential spin matrix elements: 214 ! 1 = up up = +z 215 ! 2 = down down = -z 216 ! 3 = up down = +x 217 ! 4 = down up = -y 218 Vmagconstr(ifft_local,3)=Vmagconstr(ifft_local,3) + & 219 & 2*magcon_lambda*fsm*cmm_x 220 ! & 2*magcon_lambda*fsm*(cmm_x*cos(arg)+cmm_y*sin(arg)) 221 Vmagconstr(ifft_local,4)=Vmagconstr(ifft_local,4) - & 222 & 2*magcon_lambda*fsm*cmm_y 223 ! & 2*magcon_lambda*fsm*(cmm_y*cos(arg)+cmm_x*sin(arg)) 224 ! end do loop on spirals 225 end if ! nspden 4 226 227 ! Calculate the z-component of the square bracket term 228 if (magconon==1) then 229 if (nspden == 4) then 230 ! m_z - spinat_z * <m | spinat> 231 cmm_z = intgden(4,iatom) - spinat_norm(3,iatom)*intgden_proj 232 else if (nspden == 2) then 233 ! this will be just a sign +/- : are we in the same direction as spinat_z? 234 ! need something more continuous??? To make sure the gradient pushes the state towards FM/AFM? 235 cmm_z = -sign(one, (intgden(1,iatom)-intgden(2,iatom))*spinat_norm(3,iatom)) 236 end if 237 else if (magconon==2) then 238 if (nspden == 4) then 239 cmm_z=intgden(4,iatom)-spinat(3,iatom) 240 else if (nspden == 2) then 241 ! this is up spins - down spins - requested moment ~ 0 242 ! EB: note that intgden comes from calcdensph, which, in nspden=2 case, returns 243 ! intgden(1)=rho_up=n+m 244 ! intgden(2)=rho_dn=n-m 245 ! Then, is the following line be 246 ! cmm_z=half*(intgden(1,iatom)-intgden(2,iatom)) - spinat(3,iatom) 247 ! ?? 248 cmm_z=intgden(1,iatom)-intgden(2,iatom) - spinat(3,iatom) 249 end if 250 end if 251 252 ! Calculate the constraining potential for z-component of the mag. mom. vector 253 Vmagconstr(ifft_local,1)=Vmagconstr(ifft_local,1) + 2*magcon_lambda*fsm*cmm_z 254 Vmagconstr(ifft_local,2)=Vmagconstr(ifft_local,2) - 2*magcon_lambda*fsm*cmm_z 255 ! end do loop on spirals 256 257 end do ! i1 258 end do ! i2 259 end if ! if this is my fft slice 260 end do ! i3 261 262 ! end loop over atoms 263 end do 264 265 !MPI parallelization 266 !TODO: test if xmpi_sum does the correct stuff for a slice of Vmagconstr 267 if(mpi_enreg%nproc_fft>1)then 268 call timab(48,1,tsec) 269 call xmpi_sum(Vmagconstr,mpi_enreg%comm_fft,ierr) 270 call timab(48,2,tsec) 271 end if 272 273 ! write (201,*) '# potential 1' 274 ! write (201,*) Vmagconstr(:,1) 275 276 ! write (202,*) '# potential 2' 277 ! write (202,*) Vmagconstr(:,2) 278 279 ! if (nspden > 2) then 280 ! write (203,*) '# potential 3' 281 ! write (203,*) Vmagconstr(:,3) 282 283 ! write (204,*) '# potential 4' 284 ! write (204,*) Vmagconstr(:,4) 285 ! end if 286 287 end subroutine mag_constr