what's the meaning of the tag "newg" in the code?

Total energy, geometry optimization, DFT+U, spin....

Moderator: bguster

Locked
onion2440
Posts: 34
Joined: Sat Sep 05, 2015 10:04 am

what's the meaning of the tag "newg" in the code?

Post by onion2440 » Mon Dec 14, 2015 1:56 am

Dear everyone,
I want known the meaning of the code attached below, especially the tag "newg". the code is in the abinit-7.10.5/src/41_geometry/m_ewald.F90 in about 830th row. any answers would be help for me! thank you!
yours
sincerely

newg=0
do ig3=-ng,ng
do ig2=-ng,ng
do ig1=-ng,ng
if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng .or. ng==1 )then

gpq(1)=(ig1+qphon(1))*gprim(1,1)/acell(1)+(ig2+qphon(2))*&
& gprim(1,2)/acell(2)+(ig3+qphon(3))*gprim(1,3)/acell(3)
gpq(2)=(ig1+qphon(1))*gprim(2,1)/acell(1)+(ig2+qphon(2))*&
& gprim(2,2)/acell(2)+(ig3+qphon(3))*gprim(2,3)/acell(3)
gpq(3)=(ig1+qphon(1))*gprim(3,1)/acell(1)+(ig2+qphon(2))*&
& gprim(3,2)/acell(2)+(ig3+qphon(3))*gprim(3,3)/acell(3)
gsq=0.0_dp
do jj=1,3
do ii=1,3
gpqgpq(ii,jj)=gpq(ii)*gpq(jj)
gsq=gsq+gpqgpq(ii,jj)*dielt(ii,jj)
end do
end do

! Skip q=0:
if (gsq<1.0d-20) then
if (sumg0==1) then
write(message,'(a,a,a,a,a)' )&
& 'The phonon wavelength should not be zero :',ch10,&
& 'there are non-analytical terms that cannot be treated.',ch10,&
& 'Action: subtract this wavelength from the input file.'
MSG_ERROR(message)
end if

else

arg1=(two_pi**2)*gsq/(4*eta)

! Larger arg gives 0 contribution:
if (arg1<=80._dp) then
newg=1

! Here calculate the term
term1=exp(-arg1)/gsq
do jj=1,3
do ii=1,3
gpqfac(ii,jj)=gpqgpq(ii,jj)*term1
end do
end do

! do ia=1,natom
! arga=two_pi*( (ig1+qphon(1))*xred(1,ia)&
!& +(ig2+qphon(2))*xred(2,ia)&
!& +(ig3+qphon(3))*xred(3,ia) )
! cosqxred(ia)=cos(arga)
! sinqxred(ia)=sin(arga)
! end do
! MJV: replaced old calls to cos and sin. Checked for 10 tests in v2 that max error is about 6.e-15, usually < 2.e-15
do ia=1,natom
cosqxred(ia)=real(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
sinqxred(ia)=aimag(exp2piqx(ia)*expx1(ig1, ia)*expx2(ig2, ia)*expx3(ig3, ia))
end do

! First, the diagonal terms
do nu=1,3
do ia=1,natom
do mu=nu,3
dyewt(1,mu,ia,nu,ia)=dyewt(1,mu,ia,nu,ia)+&
& gpqfac(mu,nu)
end do
end do
end do

! Then, the non-diagonal ones
do ib=2,natom
do ia=1,ib-1
c123r=cosqxred(ia)*cosqxred(ib)+sinqxred(ia)*sinqxred(ib)
c123i=sinqxred(ia)*cosqxred(ib)-cosqxred(ia)*sinqxred(ib)
! The most inner loop
do nu=1,3
do mu=nu,3
dyewt(1,mu,ia,nu,ib)=dyewt(1,mu,ia,nu,ib)+&
& gpqfac(mu,nu)*c123r
dyewt(2,mu,ia,nu,ib)=dyewt(2,mu,ia,nu,ib)+&
& gpqfac(mu,nu)*c123i
end do
end do
end do
end do

end if ! endif exp() argument is smaller than 80
end if ! Endif g/=0 :
end if ! End triple summation over Gs:
end do
end do
end do
! Check if new shell must be calculated
if(newg==0)exit
end do

Locked