TABLE OF CONTENTS
ABINIT/mkcore_inner [ Functions ]
NAME
mkcore_inner
FUNCTION
FIXME: add description.
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (T. Rangel) 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
argin(sizein)=description
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
PARENTS
mkcore_paw,mkcore_wvl
CHILDREN
paw_splint,paw_splint_der,sort_dp
SOURCE
33 #if defined HAVE_CONFIG_H 34 #include "config.h" 35 #endif 36 37 #include "abi_common.h" 38 39 subroutine mkcore_inner(corfra,core_mesh,dyfrx2,grxc1,grxc2,grxc3,ifftsph,msz,natom,ncmax,nfft,& 40 & nfgd,nfgd_r0,nspden,n3xccc,option,pawtab,rmet,rr,strdia,vxc,xccc3d,& 41 & rred) ! optional argument 42 43 use defs_basis 44 use m_profiling_abi 45 use m_errors 46 use m_sort, only : sort_dp 47 use m_pawrad, only : pawrad_type 48 use m_pawtab, only : pawtab_type 49 use m_paw_numeric, only : paw_splint,paw_splint_der 50 51 !This section has been created automatically by the script Abilint (TD). 52 !Do not modify the following lines by hand. 53 #undef ABI_FUNC 54 #define ABI_FUNC 'mkcore_inner' 55 !End of the abilint section 56 57 implicit none 58 59 !Arguments ------------------------------------ 60 !scalars 61 integer ,intent(in) :: msz,natom,ncmax,nfft,nfgd,nfgd_r0,nspden,n3xccc,option 62 real(dp),intent(out) :: grxc1,grxc2,grxc3,strdia 63 !arrays 64 integer,intent(in) :: ifftsph(ncmax) 65 real(dp),intent(in) :: rmet(3,3),rr(ncmax),vxc(nfft,nspden) 66 real(dp),intent(in),optional :: rred(:,:) 67 real(dp),intent(inout) :: corfra(3,3),dyfrx2(3,3,natom),xccc3d(n3xccc) 68 type(pawrad_type),intent(in) :: core_mesh 69 type(pawtab_type),intent(in) :: pawtab 70 71 !Local variables------------------------------- 72 !scalars 73 integer :: iatom,ifgd,ii,jj,mu,nu 74 character(len=500) :: message 75 real(dp) :: term,term2 76 !arrays 77 integer :: iindex(nfgd) 78 real(dp) :: tcore(nfgd),dtcore(nfgd),rr_tmp(nfgd) 79 real(dp),allocatable :: d2tcore(:) 80 81 ! ************************************************************************* 82 83 !Checks 84 if(nspden >1) then 85 write(message, '(a)')'mkcore_inner: this is not yet generalized to npsden>1' 86 MSG_ERROR(message) 87 end if 88 if (present(rred)) then 89 if (option>1.and.size(rred)/=3*ncmax) then 90 write(message, '(a)')'mkcore_inner: incorrect size for rred' 91 MSG_BUG(message) 92 end if 93 else if (option>1) then 94 write(message, '(a)')'mkcore_inner: rred is not present and option >1' 95 MSG_BUG(message) 96 end if 97 98 !Retrieve values of pseudo core density (and derivative) 99 rr_tmp=rr(1:nfgd) 100 iindex(1:nfgd)=(/(ii,ii=1,nfgd)/) 101 call sort_dp(nfgd,rr_tmp,iindex(1:nfgd),tol16) 102 if (option==1.or.option==3) then 103 call paw_splint(msz,core_mesh%rad,pawtab%tcoredens(:,1),pawtab%tcoredens(:,3),& 104 & nfgd,rr_tmp,tcore) 105 end if 106 if (option>=2) then 107 call paw_splint_der(msz,core_mesh%rad,pawtab%tcoredens(:,1),pawtab%tcoredens(:,3),& 108 & nfgd,rr_tmp,dtcore) 109 end if 110 111 !Accumulate contributions to core density on the entire cell 112 if (option==1) then 113 do ifgd=1,nfgd 114 ii=iindex(ifgd);jj=ifftsph(ii) 115 xccc3d(jj)=xccc3d(jj)+tcore(ifgd) 116 end do 117 118 !Accumulate contributions to Exc gradients 119 else if (option==2) then 120 do ifgd=1,nfgd 121 ii=iindex(ifgd);jj=ifftsph(ii) 122 term=vxc(jj,1)*dtcore(ifgd)/rr_tmp(ifgd) 123 grxc1=grxc1-rred(1,ii)*term 124 grxc2=grxc2-rred(2,ii)*term 125 grxc3=grxc3-rred(3,ii)*term 126 end do 127 128 !Accumulate contributions to stress tensor 129 else if (option==3) then 130 ! Write out the 6 symmetric components 131 do ifgd=1,nfgd 132 ii=iindex(ifgd);jj=ifftsph(ii) 133 term=vxc(jj,1)*dtcore(ifgd)/rr_tmp(ifgd) 134 corfra(1,1)=corfra(1,1)+term*rred(1,ifgd)**2 135 corfra(2,2)=corfra(2,2)+term*rred(2,ifgd)**2 136 corfra(3,3)=corfra(3,3)+term*rred(3,ifgd)**2 137 corfra(3,2)=corfra(3,2)+term*rred(3,ifgd)*rred(2,ifgd) 138 corfra(3,1)=corfra(3,1)+term*rred(3,ifgd)*rred(1,ifgd) 139 corfra(2,1)=corfra(2,1)+term*rred(2,ifgd)*rred(1,ifgd) 140 end do 141 ! Also compute diagonal term 142 do ifgd=1,nfgd 143 ii=iindex(ifgd);jj=ifftsph(ii) 144 strdia=strdia+vxc(jj,1)*tcore(ii) 145 end do 146 147 !Compute frozen-wf contribution to Dynamical matrix 148 else if (option==4) then 149 ABI_ALLOCATE(d2tcore,(nfgd)) 150 if(nfgd>0) then 151 ! Evaluate spline fit of 2nd der of pseudo core density 152 call paw_splint(msz,core_mesh%rad,pawtab%tcoredens(:,3),pawtab%tcoredens(:,5),& 153 & nfgd,rr_tmp,d2tcore) 154 do ifgd=1,nfgd 155 ii=iindex(ifgd);jj=ifftsph(ii) 156 term=vxc(jj,1)*dtcore(ifgd)/rr_tmp(ifgd) 157 term2=term*d2tcore(ifgd)/rr_tmp(ifgd) 158 do mu=1,3 159 do nu=1,3 160 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)& 161 & +(term2-term/rr_tmp(iindex(ifgd))**2)& 162 & *rred(mu,iatom)*rred(nu,iatom)+term*rmet(mu,nu) 163 end do 164 end do 165 end do 166 end if 167 ! Contributions from |r-R|=0 168 if (nfgd_r0>0) then 169 rr_tmp(1)=tol10 170 call paw_splint(msz,core_mesh%rad,pawtab%tcoredens(:,3),pawtab%tcoredens(:,5),& 171 & 1,rr_tmp,d2tcore(1)) 172 ifgd=1 173 ii=iindex(ifgd);jj=ifftsph(ii) 174 term=vxc(jj,1)*d2tcore(ifgd) 175 do mu=1,3 176 do nu=1,3 177 dyfrx2(mu,nu,iatom)=dyfrx2(mu,nu,iatom)+term*rmet(mu,nu) 178 end do 179 end do 180 end if 181 ABI_DEALLOCATE(d2tcore) 182 183 end if !option 184 185 end subroutine mkcore_inner