[Pw_forum] Forces and finite electric fields (via Berry phase)
Sohrab Ismail-Beigi
sohrab.ismail-beigi at yale.edu
Sat Dec 15 21:16:29 CET 2007
Dear users,
I've been experimenting with the addition of external electric fields
via the Berry phase method with pw.x --- i.e. the type of electric
field turned on by the input flag "lelfield" in conjunction with flags
"lberry, gdir, nppstr, nberrcyc, efield". After some hiccups
(regarding k-point convergence and symmetry problems) and noticing
that version 3.2.2 doesn't do this right (but the current CVS version
apparently does), it seems that most of it works reasonably with one
exception:
It seems that the atomic forces are missing the forces due to the
external field --- a simple omission but one that causes the forces to
be rather bizarre as written. I've added a short fix to "forces.f90"
inside the PW directory and I attach the modified forces.f90 file
below. But someone who knows more about the code should check that in
fact this is not being done in duplicate (although after a careful
perusing of the code, I don't see any other place the forces are added).
Has anyone else used this method successfully with pw.x and has some
comments? For example, it seems that the symmetries are handled a bit
unusually when the field is turned on --- of course, with an external
field, there are fewer symmetries but I'm not sure the ones it uses
are the "right" ones (e.g. it generally reports only one symmetry
[identity] in the output file but actually reduces the k-point set
somewhat but not as much as when there is no field turned on).
Sincerely,
---------------------------------------------
Sohrab Ismail-Beigi
Associate Professor of Applied Physics and Physics
Dept. of Applied Physics
Yale University
P.O. Box 208284
New Haven, CT 06520
Phone: 203.432.2107
Fax: 203.432.4283
Email: sohrab.ismail-beigi at yale.edu
http://volga.eng.yale.edu/ - group home page
http://www.eng.yale.edu/aphy/ - dept. page w/ links
http://www.eng.yale.edu/graduate -- info about graduate studies in AP
and brochures
P.S. The forces.f90 file is attached below. Note that this is based
on what is in the CVS file --- corrections to espresso-3.2.2 are just
the same with some minor changes to the "use" lines. My corrections
are identified by the "SIB". There is also a debugging type "write"
statement that should be removed for cleaner code...
!
! Copyright (C) 2001-2006 Quantum-ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
#include "f_defs.h"
!
!----------------------------------------------------------------------------
SUBROUTINE forces()
!----------------------------------------------------------------------------
!
! ... This routine is a driver routine which computes the forces
! ... acting on the atoms. The complete expression of the forces
! ... contains four parts which are computed by different routines:
!
! ... a) force_lc, local contribution to the forces
! ... b) force_cc, contribution due to NLCC
! ... c) force_ew, contribution due to the electrostatic ewald
term
! ... d) force_us, contribution due to the non-local potential
! ... e) force_corr, correction term for incomplete self-
consistency
! ... f) force_hub, contribution due to the Hubbard term
!
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE cell_base, ONLY : at, bg, alat, omega
USE ions_base, ONLY : nat, nsp, ityp, tau, zv, amass
USE gvect, ONLY : ngm, gstart, nr1, nr2, nr3, nrx1, nrx2,
nrx3, &
nrxx, ngl, nl, igtongl, g, gg, gcutm
USE lsda_mod, ONLY : nspin
USE symme, ONLY : s, nsym, irt
USE vlocal, ONLY : strf, vloc
USE force_mod, ONLY : force, lforce
USE scf, ONLY : rho
USE ions_base, ONLY : if_pos
USE ldaU, ONLY : lda_plus_u
USE extfield, ONLY : tefield, forcefield
USE control_flags, ONLY : gamma_only, remove_rigid_rot, lbfgs
! SIB : added support for lelfield
USE bp, ONLY : lelfield, efield, gdir
!
IMPLICIT NONE
!
REAL(DP), ALLOCATABLE :: forcenl(:,:), &
forcelc(:,:), &
forcecc(:,:), &
forceion(:,:), &
forcescc(:,:), &
forceh(:,:)
! nonlocal, local, core-correction, ewald, scf correction terms,
and hubbard
REAL(DP) :: sumfor, sumscf, bnorm, bhat(3), flelfield
INTEGER :: ipol, na
! counter on polarization
! counter on atoms
!
!
CALL start_clock( 'forces' )
!
ALLOCATE( forcenl( 3, nat ), forcelc( 3, nat ), forcecc( 3, nat ), &
forceh( 3, nat ), forceion( 3, nat ), forcescc( 3, nat ) )
!
forcescc(:,:) = 0.D0
forceh(:,:) = 0.D0
!
WRITE( stdout, '(/,5x,"Forces acting on atoms (Ry/au):", / )')
!
! ... The nonlocal contribution is computed here
!
CALL force_us( forcenl )
!
! ... The local contribution
!
CALL force_lc( nat, tau, ityp, alat, omega, ngm, ngl, igtongl, &
nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, g, rho%of_r,
nl, &
nspin, gstart, gamma_only, vloc, forcelc )
!
! ... The NLCC contribution
!
CALL force_cc( forcecc )
!
! ... The Hubbard contribution
!
IF ( lda_plus_u ) CALL force_hub( forceh )
!
! ... The ionic contribution is computed here
!
CALL force_ew( alat, nat, nsp, ityp, zv, at, bg, tau, omega, g, &
gg, ngm, gstart, gamma_only, gcutm, strf, forceion )
!
! ... The SCF contribution
!
CALL force_corr( forcescc )
!
! ... here we sum all the contributions and compute the total force
acting
! ... on the crstal
!
DO ipol = 1, 3
!
sumfor = 0.D0
!
DO na = 1, nat
!
force(ipol,na) = forcenl(ipol,na) + &
forceion(ipol,na) + &
forcelc(ipol,na) + &
forcecc(ipol,na) + &
forceh(ipol,na) + &
forcescc(ipol,na)
!
IF ( tefield ) force(ipol,na) = force(ipol,na) +
forcefield(ipol,na)
!
! SIB : if lelfield is on, add force contribution on ions
! the 2.0 factor is because e^2=(electron charge)^2=2 in
(Ryd,a0) units
IF ( lelfield ) then
bhat(:) = bg(:,gdir)
bnorm = sqrt(dot_product(bhat(:),bhat(:)))
bhat(:) = bhat(:)/bnorm
flelfield = 2.0d0*zv(ityp(na))*efield*bhat(ipol)
write( stdout , '(a,x,i1,3x,a,i5,3x,3(a,f13.9,3x))' ) &
'axis',ipol,'atom',na,'fprelel',force(ipol,na), &
'flel',flelfield,'ftot',force(ipol,na)+flelfield
force(ipol,na) = force(ipol,na) + flelfield
END IF
!
sumfor = sumfor + force(ipol,na)
!
END DO
!
! ... impose total force = 0
!
DO na = 1, nat
!
force(ipol,na) = force(ipol,na) - sumfor / DBLE( nat )
!
END DO
!
END DO
!
! ... resymmetrize (should not be needed, but ...)
!
IF ( nsym >= 1 ) THEN
!
DO na = 1, nat
CALL trnvect( force(1,na), at, bg, - 1 )
END DO
!
CALL symvect( nat, force, nsym, s, irt )
!
DO na = 1, nat
CALL trnvect( force(1,na), at, bg, 1 )
END DO
!
END IF
!
IF ( remove_rigid_rot ) &
CALL remove_tot_torque( nat, tau, amass(ityp(:)), force )
!
! ... write on output the forces
!
DO na = 1, nat
!
WRITE( stdout, 9035) na, ityp(na), force(:,na)
!
END DO
!
! ... forces on fixed coordinates are set to zero ( C.S. 15/10/2003 )
!
force(:,:) = force(:,:) * DBLE( if_pos )
forcescc(:,:) = forcescc(:,:) * DBLE( if_pos )
!
#if defined (DEBUG)
!
WRITE( stdout, '(5x,"The non-local contrib. to forces")')
DO na = 1, nat
WRITE( stdout, 9035) na, ityp(na), ( forcenl(ipol,na), ipol = 1,
3 )
END DO
WRITE( stdout, '(5x,"The ionic contribution to forces")')
DO na = 1, nat
WRITE( stdout, 9035) na, ityp(na), ( forceion(ipol,na), ipol =
1, 3 )
END DO
WRITE( stdout, '(5x,"The local contribution to forces")')
DO na = 1, nat
WRITE( stdout, 9035) na, ityp(na), ( forcelc(ipol,na), ipol = 1,
3 )
END DO
WRITE( stdout, '(5x,"The Hubbard contrib. to forces")')
DO na = 1, nat
WRITE( stdout, 9035) na, ityp(na), ( forceh(ipol,na), ipol = 1,
3 )
END DO
WRITE( stdout, '(5x,"The SCF correction term to forces")')
DO na = 1, nat
WRITE( stdout, 9035) na, ityp(na), ( forcescc(ipol,na), ipol =
1, 3 )
END DO
!
#endif
!
sumfor = 0.D0
sumscf = 0.D0
!
DO na = 1, nat
!
sumfor = sumfor + &
force(1,na)**2 + force(2,na)**2 + force(3,na)**2
!
sumscf = sumscf + &
forcescc(1,na)**2 + forcescc(2,na)**2+ forcescc(3,na)**2
!
END DO
!
sumfor = SQRT( sumfor )
sumscf = SQRT( sumscf )
WRITE( stdout, '(/5x,"Total force = ",F12.6,5X, &
& "Total SCF correction = ",F12.6)') sumfor, sumscf
!
DEALLOCATE( forcenl, forcelc, forcecc, forceh, forceion, forcescc )
!
lforce = .TRUE.
!
CALL stop_clock( 'forces' )
!
IF ( lbfgs .AND. ( sumfor < 10.D0*sumscf ) ) &
WRITE( stdout,'(5x,"SCF correction compared to forces is too large,
reduce conv_thr")')
! CALL errore( 'forces', 'scf correction on ' // &
! & 'the force is too large: reduce conv_thr', 1 )
!
RETURN
!
9035 FORMAT(5X,'atom ',I3,' type ',I2,' force = ',3F14.8)
!
END SUBROUTINE forces
More information about the users
mailing list