TABLE OF CONTENTS
ABINIT/mag_constr_e [ Functions ]
NAME
mag_constr_e
FUNCTION
This routine is called to compute the energy corresponding to constrained magnetic moments.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (ILuk) 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
magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size spinat=fixed magnetic moments vectors magcon_lambda=the size of the penalty terms
OUTPUT
Epen=penalty contribution to the total energy corresponding to the constrained potential Econstr=??? Eexp=???
PARENTS
outscfcv
CHILDREN
calcdensph,metric,wrtout
SOURCE
34 #if defined HAVE_CONFIG_H 35 #include "config.h" 36 #endif 37 38 #include "abi_common.h" 39 40 subroutine mag_constr_e(magconon,magcon_lambda,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsph,rhor,rprimd,spinat,typat,xred) 41 42 use defs_basis 43 use defs_abitypes 44 45 use m_geometry, only : metric 46 47 !This section has been created automatically by the script Abilint (TD). 48 !Do not modify the following lines by hand. 49 #undef ABI_FUNC 50 #define ABI_FUNC 'mag_constr_e' 51 use interfaces_14_hidewrite 52 use interfaces_54_abiutil 53 !End of the abilint section 54 55 implicit none 56 57 !Arguments ------------------------------------ 58 !scalars 59 integer,intent(in) :: natom,magconon,nspden,nfft,ntypat 60 real(dp),intent(in) :: magcon_lambda 61 !arrays 62 integer, intent(in) :: ngfft(18),typat(natom) 63 real(dp),intent(in) :: spinat(3,natom), rprimd(3,3) 64 real(dp),intent(in) :: ratsph(ntypat),rhor(nfft,nspden),xred(3,natom) 65 type(MPI_type),intent(in) :: mpi_enreg 66 67 !Local variables------------------------------- 68 !scalars 69 integer :: iatom,ii 70 integer :: cplex1=1 ! dummy argument for calcdensphere 71 real(dp) :: intgden_proj, Epen,Econstr,lVp, norm 72 !arrays 73 real(dp) :: intmm(3), mag_1atom(3) 74 real(dp), allocatable :: intgden(:,:) 75 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),ucvol 76 real(dp) :: spinat_norm(3,natom) 77 character(len=500) :: message 78 79 ! ********************************************************************* 80 81 !We need the metric because it is needed in calcdensph.F90 82 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 83 84 ABI_ALLOCATE (intgden, (nspden,natom)) 85 86 !We need the integrated magnetic moments 87 cplex1=1 88 call calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,std_out,ratsph,rhor,rprimd,typat,ucvol,xred,& 89 & 1,cplex1,intgden=intgden) 90 91 Epen=0 92 Econstr=0 93 lVp=0 94 95 !Loop over atoms 96 !------------------------------------------- 97 do iatom=1,natom 98 99 norm = sqrt(sum(spinat(:,iatom)**2)) 100 spinat_norm(:,iatom) = zero 101 if (norm > tol10) then 102 spinat_norm(:,iatom) = spinat(:,iatom) / norm 103 else if (magconon == 1) then 104 ! if spinat = 0 and we are imposing the direction only, skip this atom 105 cycle 106 end if 107 ! Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector 108 ! This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector 109 110 ! for the collinear spin case, set up a fictitious 3D vector along z 111 if (nspden == 4) then 112 mag_1atom(1:3) = intgden(2:4,iatom) 113 else if (nspden == 2) then 114 mag_1atom = zero 115 mag_1atom(3) = intgden(1,iatom)-intgden(2,iatom) 116 end if 117 118 intgden_proj = zero 119 intmm = zero 120 ! Calculate the square bracket term 121 if (magconon==1) then 122 intgden_proj=spinat_norm(1,iatom)*mag_1atom(1)+ & 123 & spinat_norm(2,iatom)*mag_1atom(2)+ & 124 & spinat_norm(3,iatom)*mag_1atom(3) 125 126 do ii=1,3 127 intmm(ii)=mag_1atom(ii)-spinat_norm(ii,iatom)*intgden_proj 128 end do 129 130 ! Calculate the energy Epen corresponding to the constraining potential 131 ! Econstr and lVp do not have a clear meaning (yet) 132 Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3)) 133 Econstr=Econstr-magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3)) 134 lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3)) 135 136 else if (magconon==2) then 137 do ii=1,3 138 intmm(ii)=mag_1atom(ii)-spinat(ii,iatom) 139 end do 140 141 ! Calculate the energy Epen corresponding to the constraining potential 142 ! Epen = -Econstr - lVp 143 ! Econstr = -M**2 + spinat**2 144 ! lVp = +2 M \cdot spinat 145 Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3)) 146 Econstr=Econstr-magcon_lambda*(mag_1atom(1)*mag_1atom(1)+& 147 & mag_1atom(2)*mag_1atom(2)+& 148 & mag_1atom(3)*mag_1atom(3)) & 149 & +magcon_lambda*(spinat(1,iatom)*spinat(1,iatom)+& 150 & spinat(2,iatom)*spinat(2,iatom)+& 151 & spinat(3,iatom)*spinat(3,iatom)) 152 lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3)) 153 end if 154 155 write(message, *) 'atom constraining magnetic field' 156 call wrtout(std_out,message,'COLL') 157 write(message, '(I3,A2,E12.5,A2,E12.5,A2,E12.5)') & 158 iatom,' ',magcon_lambda*intmm(1),' ',magcon_lambda*intmm(2),' ',magcon_lambda*intmm(3) 159 call wrtout(std_out,message,'COLL') 160 161 ! End loop over atoms 162 ! ------------------------------------------- 163 end do 164 165 !Printing 166 write(message, '(A17,E10.3)' ) ' magcon_lambda = ',magcon_lambda 167 call wrtout(std_out,message,'COLL') 168 write(message, '(A17,E12.5)' ) ' Lagrange penalty = ',Epen 169 call wrtout(std_out,message,'COLL') 170 write(message, '(A17,E12.5)' ) ' E_constraint = ',Econstr 171 call wrtout(std_out,message,'COLL') 172 write(message, '(A17,E12.5)' ) ' lVp = ',lVp 173 call wrtout(std_out,message,'COLL') 174 175 ABI_DEALLOCATE (intgden) 176 177 end subroutine mag_constr_e