Skip to content

ENH: 2D Polar Forces #89

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
flokno wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

ENH: 2D Polar Forces #89

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 0 additions & 11 deletions build_things.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,19 +173,8 @@ cd ../../

# now I should build a lot of codes. If we made it this far, it should be easy.
listofcodes="
dump_dynamical_matrices
phonon_dispersion_relations
crystal_structure_info
generate_structure
canonical_configuration
lineshape
samples_from_md
extract_forceconstants
atomic_distribution
pack_simulation
refine_structure
thermal_conductivity
anharmonic_free_energy
"

#some things that are not quite ok yet
Expand Down
97 changes: 97 additions & 0 deletions log.txt
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion src/extract_forceconstants/options.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,7 @@ subroutine parse(opts)
if (lo_status .ne. 0) stop
call cli%add(switch='--polarcorrectiontype', switch_ab='-pc', hidden=.true., &
help='What kind of polar correction to use.', &
required=.false., act='store', def='3', choices='1,2,3', error=lo_status)
required=.false., act='store', def='3', choices='1,2,3,4', error=lo_status)
if (lo_status .ne. 0) stop
call cli%add(switch='--secondorder_njump', switch_ab='-nj2', &
help='Second order neighbour jumps.', &
Expand Down
9 changes: 7 additions & 2 deletions src/libolle/ifc_solvers_diel.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -498,8 +498,13 @@ module subroutine lo_solve_for_borncharges(map, p, Z, eps, filename, mw, mem, ve
! Make sure it's hermitian? Seems like the sensible thing to do.
call lo_coeffmatrix_unitcell_Z_singlet(map, coeffM)

call ew%set(p, epsbuf, 2, 1E-20_r8, verbosity)
call ew%force_borncharges_Hermitian(map%xuc%x_Z_singlet, coeffM, epsbuf, p, verbosity)
if (map%polarcorrectiontype .eq. 4) then
call ew%set_2D(p, epsbuf, -1, 1E-20_r8, verbosity)
call ew%force_borncharges_Hermitian(map%xuc%x_Z_singlet, coeffM, epsbuf, p, verbosity)
else
call ew%set_2D(p, epsbuf, -1, 1E-20_r8, verbosity)
call ew%force_borncharges_Hermitian(map%xuc%x_Z_singlet, coeffM, epsbuf, p, verbosity)
endif

! cleanup
call mem%deallocate(coeffM, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
Expand Down
36 changes: 33 additions & 3 deletions src/libolle/ifc_solvers_prepsolver.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -437,11 +437,18 @@ module subroutine coefficients_ifc(ih, tp, map, sim, uc, ss, mw, mem, verbosity)
end if

! If it's polar, it might be a good idea to be able to remove those forces.
if (map%polar .gt. 0 .and. map%polarcorrectiontype .eq. 3) then
if (map%polar .eq. 1) then
! Get some forceconstant thingy. This is a little confusing since
! things can be defined in different ways.
call map%get_secondorder_forceconstant(uc, fc, mem, verbosity=-1)
if (fc%polar) then

if (fc%polar .and. map%polarcorrectiontype .eq. 3) then
allocate (ih%polar_fc(3, 3, ss%na, ss%na))
ih%polar_fc = 0.0_r8
call fc%supercell_longrange_dynamical_matrix_at_gamma(ss, ih%polar_fc, 1E-15_r8)
else if (fc%polar .and. map%polarcorrectiontype .eq. 4) then
! AA: polar forces for 2D materials
write(*,*) ('Getting the LR IFCS in 2D to subtract the forces: ')
allocate (ih%polar_fc(3, 3, ss%na, ss%na))
ih%polar_fc = 0.0_r8
call fc%supercell_longrange_dynamical_matrix_at_gamma(ss, ih%polar_fc, 1E-15_r8)
Expand All @@ -452,7 +459,11 @@ module subroutine coefficients_ifc(ih, tp, map, sim, uc, ss, mw, mem, verbosity)
else
allocate (ih%polar_fc(1, 1, 1, 1))
ih%polar_fc = -lo_huge

end if



end block init

coeff: block
Expand Down Expand Up @@ -503,7 +514,8 @@ module subroutine coefficients_ifc(ih, tp, map, sim, uc, ss, mw, mem, verbosity)
end if

! And perhaps the polar force thingy
if (map%polar .gt. 0 .and. map%polarcorrectiontype .eq. 3 .and. map%xuc%nx_Z_singlet .gt. 0) then
! AA: do not force the correction to be type 3!
if (map%polar .gt. 0 .and. map%xuc%nx_Z_singlet .gt. 0) then
pf = 0.0_r8
epol = 0.0_r8
do j = 1, ss%na
Expand All @@ -516,6 +528,23 @@ module subroutine coefficients_ifc(ih, tp, map, sim, uc, ss, mw, mem, verbosity)
pf = 0.0_r8
epol = 0.0_r8
end if

!write(*,*) ('I want to write the polar forces out:')

!if (tt == 1) then
! open(unit=10, file='pf', status='replace', action='write')
!else
! open(unit=10, file='pf', status='old', position='append', action='write')
!endif

!do j = 1, ss%na
! write(10, '(3F10.5)') pf(:, j)
!end do

!close(10)




! Store
ii = (tt - 1)*nf + 1
Expand All @@ -542,6 +571,7 @@ module subroutine coefficients_ifc(ih, tp, map, sim, uc, ss, mw, mem, verbosity)
end if
end do


if (map%have_fc_singlet) call mem%deallocate(CMD1, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
if (map%have_fc_pair) call mem%deallocate(CMD2, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
if (map%have_fc_triplet) call mem%deallocate(CMD3, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__)
Expand Down
183 changes: 182 additions & 1 deletion src/libolle/lo_longrange_electrostatics.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,12 +34,18 @@ module lo_longrange_electrostatics
contains
!> set parameters
procedure :: set => lo_set_dipole_ewald_parameters
!> set parameters for 2D
procedure :: set_2D => lo_set_dipole_ewald_parameters_2D
!> force Born charges Hermitian
procedure :: force_borncharges_hermitian
!> calculate the long-range dynamical matrix
procedure :: longrange_dynamical_matrix
!> calculate the long-range 2D dynamical matrix
procedure :: longrange_dynamical_matrix_2D
!> calculate the long-range forceconstant for a supercell
procedure :: supercell_longrange_forceconstant
!> calculate the long-range forceconstant for a 2D supercell
procedure :: supercell_longrange_forceconstant_2D
!> measure size in memoryx
procedure :: size_in_mem
!> destroy
Expand Down Expand Up @@ -68,6 +74,19 @@ module subroutine longrange_dynamical_matrix( &
logical, intent(in), optional :: reconly
logical, intent(in), optional :: chgmult
end subroutine
module subroutine longrange_dynamical_matrix_2D( &
ew, p, q, born_effective_charges, born_onsite_correction, eps, D, Dx, Dy, Dz, reconly, chgmult)
class(lo_ewald_parameters), intent(in) :: ew
type(lo_crystalstructure), intent(in) :: p
real(r8), dimension(3), intent(in) :: q
real(r8), dimension(:, :, :), intent(in) :: born_effective_charges
real(r8), dimension(:, :, :), intent(in) :: born_onsite_correction
real(r8), dimension(3, 3), intent(in) :: eps
complex(r8), dimension(:, :, :, :), intent(out) :: D
complex(r8), dimension(:, :, :, :), intent(out), optional :: Dx, Dy, Dz
logical, intent(in), optional :: reconly
logical, intent(in), optional :: chgmult
end subroutine
module subroutine supercell_longrange_forceconstant(ew, born_effective_charges, eps, ss, forceconstant, thres)
class(lo_ewald_parameters), intent(in) :: ew
real(r8), dimension(:, :, :), intent(in) :: born_effective_charges
Expand All @@ -76,6 +95,14 @@ module subroutine supercell_longrange_forceconstant(ew, born_effective_charges,
real(r8), dimension(:, :, :, :), intent(out) :: forceconstant
real(r8), intent(in) :: thres
end subroutine
module subroutine supercell_longrange_forceconstant_2D(ew, born_effective_charges, eps, ss, forceconstant, thres)
class(lo_ewald_parameters), intent(in) :: ew
real(r8), dimension(:, :, :), intent(in) :: born_effective_charges
real(r8), dimension(3, 3), intent(in) :: eps
type(lo_crystalstructure), intent(in) :: ss
real(r8), dimension(:, :, :, :), intent(out) :: forceconstant
real(r8), intent(in) :: thres
end subroutine
module pure function ewald_H_thingy(x, y, inveps) result(H)
real(r8), dimension(3), intent(in) :: x
real(r8), intent(in) :: y
Expand Down Expand Up @@ -194,7 +221,9 @@ subroutine lo_set_dipole_ewald_parameters(ew, p, eps, strategy, tol, verbosity,

! Decide on a sensible distance where things should be zero, I heuristically
! pick double the nearest neighbour distance. Not a crucial choice, but it works.
doublenndist = p%nearest_neighbour_distance()*2.0_r8
doublenndist = p%nearest_neighbour_distance()*2.0_r8
!doublenndist = 4.60_r8
!write(*,*) ('This is the distance where things are zero: '),doublenndist
! inverse of the metric
invdete = 1.0_r8/sqrt(lo_determ(eps))

Expand Down Expand Up @@ -385,6 +414,158 @@ subroutine lo_set_dipole_ewald_parameters(ew, p, eps, strategy, tol, verbosity,
end if
end subroutine

!> pick optimum Ewald parameters for dipole terms
subroutine lo_set_dipole_ewald_parameters_2D(ew, p, eps, strategy, tol, verbosity, forced_lambda)
!> ewald parameters
class(lo_ewald_parameters), intent(out) :: ew
!> crystal structure
type(lo_crystalstructure), intent(in) :: p
!> dielectric tensor
real(r8), dimension(3, 3), intent(in) :: eps
!> strategy to choose lambda
integer, intent(in) :: strategy
!> tolerance
real(r8), intent(in) :: tol
!> talk a lot
integer, intent(in) :: verbosity
!> force a certain lambda parameter
real(r8), intent(in), optional :: forced_lambda

integer, parameter :: npts = 40
real(r8), dimension(3, npts) :: pts
real(r8) :: t0, t1, L, pol_out

t0 = walltime()
t1 = t0

! Set some simple things
init: block
ew%eps = eps
end block init

!AA: for now I dont optimize L, i think if we want to do that we need to rewire the code again

L = norm2(p%latticevectors(:,3))
pol_out = (L/(4.0_r8*lo_pi))*(1 - 1/eps(3,3)) !out-of-plane polarizability

ew%lambda = 1.20_r8*4*lo_pi*pol_out !setting (L) slightly higher than the stability criteria!

!ew%lambda = 4.50_r8 !setting (L) slightly higher than the stability criteria!




if (verbosity .gt. 0) then
t1 = walltime()
write (lo_iou, *) '... decided on seperation parameter (L)=', tochar(ew%lambda), ' (', tochar(t1 - t0), 's)'
t0 = t1
end if

! Now build the vectors
buildvecs: block
integer, parameter :: niter = 1000
integer, dimension(3) :: bd
real(r8), dimension(:, :), allocatable :: dr
real(r8), dimension(:), allocatable :: vn
real(r8), dimension(3, 3) :: m0
real(r8), dimension(3) :: v0
real(r8) :: krad, rrad, f0, f1
integer, dimension(:), allocatable :: di
integer :: i, j, k, l, ndim

krad = (1/ew%lambda)*10.0_r8 !the inverse of the seperation parameter
!krad = ew%lambda !the inverse of the seperation parameter

! Then exactly the same thing again for G-vectors
!ndim = 0
!AA: assume this is fine for now
!do
! ndim = ndim + 1
! do i = 1, 3
! m0(:, i) = p%reciprocal_latticevectors(:, i)*(2*ndim + 1)
! end do
! f0 = lo_inscribed_sphere_in_box(m0)
! if (f0 .gt. krad) exit

!end do

bd = 1
do
do i = 1, 3
m0(:, i) = p%reciprocal_latticevectors(:, i)*(2*bd(i) + 1)
end do
f0 = lo_inscribed_sphere_in_box(m0)
if (f0 .gt. krad) exit
bd = increment_dimensions(bd, p%reciprocal_latticevectors)
end do



write(*,*) ('in the ewald module, here is the supercell'),m0(:,:)
write(*,*) ('in the ewald module, here is the ndim'),bd(:)
write(*,*) ('in the ewald module, here is the krad'),krad


! count realspace vectors
f0 = krad**2
l = 0
do i = -bd(1), bd(1)
do j = -bd(2), bd(2)
v0 = p%fractional_to_cartesian(real([i, j, 0], r8), reciprocal=.true.)
if (lo_sqnorm(v0) .lt. f0) then
l = l + 1
end if
end do
end do
ew%n_Gvector = l
allocate (ew%Gvec(3, l))
allocate (dr(3, l))
allocate (vn(l))
allocate (di(l))
ew%Gvec = 0.0_r8
dr = 0.0_r8
vn = 0.0_r8
di = 0
l = 0
do i = -bd(1), bd(1)
do j = -bd(2), bd(2)
v0 = p%fractional_to_cartesian(real([i, j, 0], r8), reciprocal=.true.)
f1 = lo_sqnorm(v0)
if (f1 .lt. f0) then
l = l + 1
vn(l) = -f1
dr(:, l) = v0
end if
end do
end do
! Sort backwards by length
call lo_qsort(vn, di)
do i = 1, ew%n_Gvector
ew%Gvec(:, i) = dr(:, di(i))
end do
deallocate (dr)
deallocate (vn)
deallocate (di)

write(*,*) ('In the ewald module, Getting the Gvecs in 2D: ')
open(unit=10, file='Gvecs_ewald', status='replace', action='write')
do i = 1, ew%n_Gvector
write(10, '(3F10.5)') ew%Gvec(:, i)*lo_twopi
end do
close(10)
write(*,*) ('In the Ewald module , here is n_gvecs'),ew%n_Gvector

end block buildvecs


if (verbosity .gt. 0) then
t1 = walltime()
write (lo_iou, *) '... stored ', tochar(ew%n_Gvector), &
' reciprocal lattice vectors (', tochar(t1 - t0), 's)'
t0 = t1
end if
end subroutine

!> turn the ratio of k-space volume to r-space volume into a function to be minimized
subroutine ewald_dipole_k_r(lambda, pts, eps, inveps, p, tol, rrad, krad)
!> ewald parameter
Expand Down
Loading