Octopus
stress.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2016 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!! Copyright (C) 2023 N. Tancogne-Dejean
3!!
4!! This program is free software; you can redistribute it and/or modify
5!! it under the terms of the GNU General Public License as published by
6!! the Free Software Foundation; either version 2, or (at your option)
7!! any later version.
8!!
9!! This program is distributed in the hope that it will be useful,
10!! but WITHOUT ANY WARRANTY; without even the implied warranty of
11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12!! GNU General Public License for more details.
13!!
14!! You should have received a copy of the GNU General Public License
15!! along with this program; if not, write to the Free Software
16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17!! 02110-1301, USA.
18!!
19
20#include "global.h"
21
22! ---------------------------------------------------------
25module stress_oct_m
28 use comm_oct_m
29 use debug_oct_m
33 use energy_oct_m
35 use epot_oct_m
36 use global_oct_m
37 use grid_oct_m
39 use io_oct_m
41 use ions_oct_m
43 use, intrinsic :: iso_fortran_env
47 use lda_u_oct_m
50 use math_oct_m
51 use mesh_oct_m
55 use mpi_oct_m
60 use ps_oct_m
62 use space_oct_m
72 use types_oct_m
73 use unit_oct_m
75 use v_ks_oct_m
77 use xc_oct_m
78 use xc_f03_lib_m
80 implicit none
81
82 private
83 public :: &
87
88contains
89
90 ! ---------------------------------------------------------
92 subroutine stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
93 type(namespace_t), intent(in) :: namespace
94 type(grid_t), intent(inout) :: gr
95 type(hamiltonian_elec_t), intent(inout) :: hm
96 type(states_elec_t), target, intent(inout) :: st
97 type(ions_t), intent(inout) :: ions
98 type(v_ks_t), intent(in) :: ks
99 type(partner_list_t), intent(in) :: ext_partners
100
101 real(real64), allocatable :: rho_total(:)
102 real(real64) :: stress(3,3) ! stress tensor in Cartesian coordinate
103 real(real64) :: stress_kin(3,3), stress_Hartree(3,3), stress_xc(3,3), stress_xc_nlcc(3,3)
104 real(real64) :: stress_ps(3,3), stress_ps_nl(3,3), stress_ps_local(3,3), stress_ii(3,3)
105 real(real64) :: stress_hubbard(3,3)
106 integer :: ip
107 real(real64), allocatable :: vh(:)
108 real(real64), allocatable :: grad_vh(:,:)
109 real(real64) :: ehartree
110 real(real64), contiguous, pointer :: rho(:)
111
112 call profiling_in("STRESS_CALCULATE")
113 push_sub(stress_calculate)
114
115 if (st%wfs_type /= type_cmplx) then
116 write(message(1),'(a)') 'The stress tensors for real wavefunctions has not been implemented!'
117
118 if (hm%kpoints%full%npoints == 1) then
119 write(message(2),'(a)') 'For testing this feature, you can add ForceComplex=yes to the input file'
120 call messages_fatal(2, namespace=namespace)
121 end if
122
123 call messages_fatal(1, namespace=namespace)
124 end if
125
126 if (ions%space%periodic_dim == 1) then
127 call messages_not_implemented("Stress tensor for 1D periodic systems", namespace=namespace)
128 end if
129
130 if (.not. ions%space%is_periodic()) then
131 write(message(1),'(a)') 'The stress tensor cannot be computed for isolated systems'
132 call messages_fatal(1, namespace=namespace)
133 end if
134
135 if (ks%vdw%vdw_correction /= option__vdwcorrection__none .and. .not. any(ks%vdw%vdw_correction == d3_lib_options)) then
136 write(message(1),'(a)') 'The stress tensor is currently only implemented with DFT-D3 vdW correction'
137 call messages_fatal(1, namespace=namespace)
138 end if
139
140 if (hm%pcm%run_pcm) then
141 call messages_not_implemented('Stress tensor with PCM')
142 end if
143
144 if (allocated(hm%v_static)) then
145 call messages_not_implemented('Stress tensor with static electric fields')
146 end if
147
148 if (ks%has_photons) then
149 call messages_not_implemented('Stress tensor with photon modes')
150 end if
151
152 if (.not. hm%vnl%apply_projector_matrices) then
153 call messages_not_implemented('Stress tensor with relativistic Kleinman-Bylander pseudopotential')
154 end if
155
156 if (hm%ep%reltype == scalar_relativistic_zora .or. hm%ep%reltype == fully_relativistic_zora) then
157 call messages_not_implemented('Stress tensor with ZORA')
158 end if
159
160 ! Checks for the xc part of KS-DFT and GKS-DFT
161 if (ks%theory_level == kohn_sham_dft .or. ks%theory_level == generalized_kohn_sham_dft) then
162 if (.not. xc_is_energy_functional(hm%xc)) then
163 call messages_not_implemented("Stress tensor with xc functionals that are not energy functionals")
164 end if
165
166 if ( .not. in_family(hm%xc%family, [xc_family_lda, xc_family_gga])) then
167 write(message(1),'(a)') 'The stress tensor computation is currently only possible at the Kohn-Sham DFT level'
168 write(message(2),'(a)') 'with LDA and GGA functionals or for independent particles.'
169 call messages_fatal(2, namespace=namespace)
170 end if
171
172 if (in_family(hm%xc%family, [xc_family_gga]) .and. st%d%ispin == spinors) then
173 call messages_not_implemented("Stress tensor for GGAs with spinors", namespace=namespace)
174 end if
175 end if
176
177 if (hm%magnetic_constrain%level /= constrain_none) then
178 call messages_not_implemented("Stress tensor with MagneticConstrain /= constrain_none")
179 end if
180
181 stress(:,:) = m_zero
182
183 safe_allocate(rho_total(1:gr%np_part))
184 do ip = 1, gr%np
185 rho_total(ip) = sum(st%rho(ip, 1:st%d%nspin))
186 end do
188 ! As we rely on some of the full energy components, we need to recompute it first
189 ! TODO: We should restrict the components of the energy needed to be computed
190 call energy_calc_total(namespace, ions%space, hm, gr, st, ext_partners, iunit = -1, full = .true.)
191
192 ! In order to get the electrostatic part (Hartree and local pseudopotential part),
193 ! we need to get the Hartree potential and its gradient
194 safe_allocate(vh(1:gr%np_part))
195 safe_allocate(grad_vh(1:gr%np, 1:gr%der%dim))
196 if (ks%theory_level /= independent_particles) then
197 call lalg_copy(gr%np, hm%ks_pot%vhartree, vh)
198 else
199 if (hm%d%spin_channels > 1) then
200 safe_allocate(rho(1:gr%np_part))
201 call lalg_copy(gr%np, st%rho(:,1), rho)
202 call lalg_axpy(gr%np, m_one, st%rho(:,2), rho)
203 else
204 rho => st%rho(:,1)
205 end if
206 ! In the case of independent particles, we use the electron density without NLCC
207 call dpoisson_solve(hm%psolver, ions%namespace, vh, rho, all_nodes = .true.)
208 if (hm%d%spin_channels > 1) then
209 safe_deallocate_p(rho)
210 else
211 nullify(rho)
212 end if
213 end if
214 ehartree = hm%energy%hartree
215 ! We also compute the gradient here
216 call dderivatives_grad(gr%der, vh, grad_vh)
217
218 ! We now compute the various contributions to the stress tensor
219
220 ! Stress from kinetic energy of electrons
221 call stress_from_kinetic(gr, ions%space, hm, st, gr%symm, ions%latt%rcell_volume, stress_kin)
222 stress = stress + stress_kin
223
224 if (ks%theory_level == independent_particles) then
225 stress_hartree = m_zero
226 stress_xc = m_zero
227 else
228 call stress_from_hartree(gr, ions%space, ions%latt%rcell_volume, grad_vh, ehartree, stress_hartree)
229 stress = stress + stress_hartree
230
231 call stress_from_xc(hm%energy, ions%latt%rcell_volume, ions%space%periodic_dim, stress_xc)
232
233 ! Nonlinear core correction contribution
234 if (allocated(st%rho_core)) then
235 call stress_from_xc_nlcc(ions%latt%rcell_volume, gr, st, ions, hm%ks_pot%vxc, stress_xc_nlcc)
236 stress_xc = stress_xc + stress_xc_nlcc
237 end if
238 ! Adds the beyond LDA contribution to the stress tensor
239 stress_xc = stress_xc + ks%stress_xc_gga / ions%latt%rcell_volume
240 stress = stress + stress_xc
241 end if
242
243 call stress_from_pseudo_local(gr, hm, ions, rho_total, grad_vh, stress_ps_local)
244 stress_ps = stress_ps_local
245 stress = stress + stress_ps_local
246
247 safe_deallocate_a(vh)
248 safe_deallocate_a(grad_vh)
249
250 call stress_from_pseudo_nonloc(gr, st, hm, ions, stress_ps_nl)
251 stress_ps = stress_ps + stress_ps_nl
252 stress = stress + stress_ps_nl
253
254 call stress_from_hubbard(namespace, gr, st, hm, ions%space, ions%latt%rcell_volume, stress_hubbard)
255 stress = stress + stress_hubbard
256
257 call ion_interaction_stress(ions%ion_interaction, ions%space, ions%latt, ions%atom, ions%natoms, ions%pos, stress_ii)
258 stress = stress + stress_ii
259 ! Stress from kinetic energy of ion
260 ! Stress from ion-field interaction
261
262 ! Sign changed to fit conventional definition
263 stress = -stress
264
265 st%stress_tensors%kinetic = stress_kin
266 st%stress_tensors%Hartree = stress_hartree
267 st%stress_tensors%xc = stress_xc
268 st%stress_tensors%ps_local = stress_ps_local
269 st%stress_tensors%ps_nl = stress_ps_nl
270 st%stress_tensors%hubbard = stress_hubbard
271 st%stress_tensors%ion_ion = stress_ii
272
273 ! Stress contribution from vdW D3
274 if (ks%vdw%vdw_correction /= option__vdwcorrection__none) then
275 st%stress_tensors%vdw = hm%ep%vdw_stress
276 else
277 st%stress_tensors%vdw = m_zero
278 end if
279 stress = stress + st%stress_tensors%vdw
280
281 ! Symmetrize the stress tensor if we use k-point symmetries
282 if (hm%kpoints%use_symmetries) then
283 call dsymmetrize_tensor_cart(gr%symm, stress, use_non_symmorphic=.true.)
284 end if
285 ! We guarantee that the matrix is truely symmetric. There could be small numerical assymetries after symmetrization
286 call symmetrize_matrix(ions%space%periodic_dim, stress)
287
288 st%stress_tensors%total = stress
289
290 ! Some sumrule for validation
291 ! Sumrule is -3P_{kin}\Omega = 2 E_{kin}
292 st%stress_tensors%kinetic_sumrule = m_zero
293 ! Sumrule is -3P_{Hartree}\Omega = E_{Hartree}
294 st%stress_tensors%Hartree_sumrule = m_zero
295 if(ions%space%periodic_dim == 3) then
296 st%stress_tensors%kinetic_sumrule = (stress_kin(1,1) + stress_kin(2,2) + stress_kin(3,3))*ions%latt%rcell_volume
297 st%stress_tensors%kinetic_sumrule = st%stress_tensors%kinetic_sumrule - m_two * hm%energy%kinetic
298
299 st%stress_tensors%hartree_sumrule = (stress_hartree(1,1) + stress_hartree(2,2) + stress_hartree(3,3))*ions%latt%rcell_volume
300 st%stress_tensors%hartree_sumrule = st%stress_tensors%hartree_sumrule - hm%energy%hartree
301 end if
302
303 safe_deallocate_a(rho_total)
304
305 pop_sub(stress_calculate)
306 call profiling_out("STRESS_CALCULATE")
307 end subroutine stress_calculate
308
309 ! -------------------------------------------------------
324 subroutine stress_from_kinetic(gr, space, hm, st, symm, rcell_volume, stress_kin)
325 type(grid_t), intent(in) :: gr
326 class(space_t), intent(in) :: space
327 type(hamiltonian_elec_t), intent(in) :: hm
328 type(states_elec_t), intent(inout) :: st
329 type(symmetries_t), intent(in) :: symm
330 real(real64), intent(in) :: rcell_volume
331 real(real64), intent(out) :: stress_kin(3, 3)
332
333 integer :: ik, ist, idir, jdir, ib, minst, maxst
334 complex(real64), allocatable :: stress_l_block(:)
335 type(wfs_elec_t) :: psib, gpsib(space%dim)
336
337 call profiling_in("STRESS_FROM_KINETIC")
338 push_sub(stress_from_kinetic)
339
340 stress_kin(:,:) = m_zero
341
342 safe_allocate(stress_l_block(1:st%block_size))
343
344 do ik = st%d%kpt%start, st%d%kpt%end
345 if (st%kweights(ik) <= m_epsilon) cycle
346
347 do ib = st%group%block_start, st%group%block_end
348 minst = states_elec_block_min(st, ib)
349 maxst = states_elec_block_max(st, ib)
350
351 call hamiltonian_elec_copy_and_set_phase(hm, gr, st%d%kpt, st%group%psib(ib, ik), psib)
352
353 ! calculate the gradient
354 call zderivatives_batch_grad(gr%der, psib, gpsib, set_bc=.false.)
355
356 ! Accumulate the result
357 do idir = 1, space%periodic_dim
358 do jdir = idir, space%periodic_dim
359 call zmesh_batch_dotp_vector(gr, gpsib(idir), gpsib(jdir), stress_l_block)
360
361 do ist = minst, maxst
362 stress_kin(idir,jdir) = stress_kin(idir,jdir) &
363 + st%kweights(ik) * st%occ(ist, ik) &
364 * real(stress_l_block(ist - minst + 1), real64)
365 end do
366 end do
367 end do
368
369 do idir = 1, space%dim
370 call gpsib(idir)%end()
371 end do
372 call psib%end()
373
374 end do
375 end do
376
377 if (st%parallel_in_states .or. st%d%kpt%parallel) then
378 call comm_allreduce(st%st_kpt_mpi_grp, stress_kin)
379 end if
380
381
382 ! Symmetrize the kinetic stress tensor
383 call upper_triangular_to_hermitian(space%periodic_dim, stress_kin)
384
385 ! Symmetrize the stress tensor if we use k-point symmetries
386 if (hm%kpoints%use_symmetries) then
387 call dsymmetrize_tensor_cart(symm, stress_kin, use_non_symmorphic=.true.)
388 end if
389
390 stress_kin = stress_kin / rcell_volume
391
392 call profiling_out("STRESS_FROM_KINETIC")
393 pop_sub(stress_from_kinetic)
394 end subroutine stress_from_kinetic
395
396 ! -------------------------------------------------------
413 subroutine stress_from_hartree(gr, space, volume, grad_vh, ehartree, stress_Hartree)
414 type(grid_t), intent(in) :: gr
415 class(space_t), intent(in) :: space
416 real(real64), intent(in) :: volume
417 real(real64), intent(in) :: grad_vh(:,:)
418 real(real64), intent(in) :: ehartree
419 real(real64), intent(out) :: stress_Hartree(3, 3)
420
421 integer :: idir, jdir
422
423 call profiling_in("STRESS_FROM_HARTREE")
424 push_sub(stress_from_hartree)
425
426 stress_hartree(:,:) = m_zero
427
428 do idir = 1, space%periodic_dim
429 do jdir = idir, space%periodic_dim
430 stress_hartree(idir, jdir) = -dmf_dotp(gr, grad_vh(:,idir), grad_vh(:, jdir))/m_four/m_pi
431 end do
432 stress_hartree(idir, idir) = stress_hartree(idir, idir) + ehartree
433 end do
434
435 call upper_triangular_to_hermitian(space%periodic_dim, stress_hartree)
436
437 stress_hartree = stress_hartree/volume
438
439 call profiling_out("STRESS_FROM_HARTREE")
440 pop_sub(stress_from_hartree)
441 end subroutine stress_from_hartree
442
443
444 ! -------------------------------------------------------
458 !
459 ! Note: We assume hm%energy%echange, correlation, and intnvxc
460 ! have already been calculated somewhere else.
461 subroutine stress_from_xc(energy, rcell_volume, periodic_dim, stress_xc)
462 type(energy_t), intent(in) :: energy
463 real(real64), intent(in) :: rcell_volume
464 integer, intent(in) :: periodic_dim
465 real(real64), intent(out) :: stress_xc(3, 3)
466
467 integer :: idir
468
469 call profiling_in("STRESS_FROM_XC")
470 push_sub(stress_from_xc)
471
472 stress_xc = m_zero
473 do idir = 1, periodic_dim
474 stress_xc(idir, idir) = - energy%exchange - energy%correlation + energy%intnvxc
475 end do
476 stress_xc(:,:) = stress_xc(:,:) / rcell_volume
477
478 call profiling_out("STRESS_FROM_XC")
479 pop_sub(stress_from_xc)
480 end subroutine stress_from_xc
481
482
483 ! -------------------------------------------------------
491 subroutine stress_from_xc_nlcc(rcell_volume, gr, st, ions, vxc, stress_xc_nlcc)
492 real(real64), intent(in) :: rcell_volume
493 type(grid_t), intent(in) :: gr
494 type(states_elec_t), intent(in) :: st
495 type(ions_t), intent(in) :: ions
496 real(real64), intent(in) :: vxc(:,:)
497 real(real64), intent(out) :: stress_xc_nlcc(3, 3)
498
499 integer :: idir, jdir, iat
500 real(real64), allocatable :: gnlcc(:,:), gnlcc_x(:,:,:), vxc_tot(:)
501
502 call profiling_in("STRESS_FROM_XC_NLCC")
503 push_sub(stress_from_xc_nlcc)
504
505 assert(allocated(st%rho_core))
506
507 stress_xc_nlcc = m_zero
509 ! We first accumulate the contribution from all the pseudo-ions
510 safe_allocate(gnlcc(gr%np, gr%der%dim))
511 safe_allocate(gnlcc_x(gr%np, gr%der%dim, gr%der%dim))
512 gnlcc_x = m_zero
513 do iat = ions%atoms_dist%start, ions%atoms_dist%end
514 assert(ions%atom(iat)%species%is_ps())
515 call species_get_nlcc_grad(ions%atom(iat)%species, ions%space, ions%latt, &
516 ions%pos(:,iat), gr, gnlcc, gnlcc_x)
517 end do
518 safe_deallocate_a(gnlcc)
519
520 if (ions%atoms_dist%parallel) then
521 call comm_allreduce(ions%atoms_dist%mpi_grp, gnlcc_x)
522 end if
523
524 ! Sum over spin of the xc potential
525 safe_allocate(vxc_tot(1:gr%np))
526 call lalg_copy(gr%np, vxc(:, 1), vxc_tot)
527 if(st%d%nspin > 1) call lalg_axpy(gr%np, m_one, vxc(:, 2), vxc_tot)
528
529 do idir = 1, ions%space%periodic_dim
530 do jdir = idir, ions%space%periodic_dim
531 stress_xc_nlcc(idir, jdir) = dmf_dotp(gr, vxc_tot, gnlcc_x(:,idir, jdir))
532 end do
533 end do
534 safe_deallocate_a(vxc_tot)
535 safe_deallocate_a(gnlcc_x)
536
537 call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_xc_nlcc)
538
539 stress_xc_nlcc(:,:) = stress_xc_nlcc(:,:) / rcell_volume
540
541 call profiling_out("STRESS_FROM_XC_NLCC")
542 pop_sub(stress_from_xc_nlcc)
543 end subroutine stress_from_xc_nlcc
544
545 ! -------------------------------------------------------
565 subroutine stress_from_pseudo_nonloc(gr, st, hm, ions, stress_ps_nl)
566 type(grid_t), target, intent(in) :: gr
567 type(states_elec_t), intent(inout) :: st
568 type(hamiltonian_elec_t), intent(in) :: hm
569 type(ions_t), intent(in) :: ions
570 real(real64), intent(out) :: stress_ps_nl(3, 3)
571
572 integer :: ik, ist, idir, jdir
573 integer :: ib, minst, maxst
574 type(wfs_elec_t) :: psib, rvnl_psib(3), gpsib(3)
575 complex(real64), allocatable :: stress_tmp(:)
576
577 call profiling_in("STRESS_FROM_PSEUDO_NL")
579
580 assert(st%wfs_type == type_cmplx)
581
582 safe_allocate(stress_tmp(1:st%block_size))
583
584 stress_ps_nl = m_zero
585
586 do ik = st%d%kpt%start, st%d%kpt%end
587
588 if (st%kweights(ik) <= m_epsilon) cycle
589
590 do ib = st%group%block_start, st%group%block_end
591 minst = states_elec_block_min(st, ib)
592 maxst = states_elec_block_max(st, ib)
593
594 call hamiltonian_elec_copy_and_set_phase(hm, gr, st%d%kpt, st%group%psib(ib, ik), psib)
595
596 ! calculate the gradient
597 call zderivatives_batch_grad(gr%der, psib, gpsib, set_bc=.false.)
598
599
600 ! Get rV_NL |\psi> for all atoms
601 do idir = 1, gr%der%dim
602 call psib%copy_to(rvnl_psib(idir))
603 call batch_set_zero(rvnl_psib(idir))
604 end do
605 call hm%vnl%zr_vn_local(gr, st%d, gr%der%boundaries%spiral, psib, rvnl_psib)
606
607 do idir = 1, ions%space%periodic_dim
608 do jdir = idir, ions%space%periodic_dim
609 call zmesh_batch_dotp_vector(gr, gpsib(idir), rvnl_psib(jdir), stress_tmp)
610
611 do ist = minst, maxst
612 stress_ps_nl(idir, jdir) = stress_ps_nl(idir, jdir) &
613 + m_two * st%kweights(ik) * st%occ(ist, ik) * real(stress_tmp(ist-minst+1), real64)
614 end do
615
616 end do
617 end do
618
619 do idir = 1, gr%der%dim
620 call rvnl_psib(idir)%end()
621 call gpsib(idir)%end()
622 end do
623 call psib%end()
624 end do
625 end do
626
627 safe_deallocate_a(stress_tmp)
628
629 if (st%parallel_in_states .or. st%d%kpt%parallel) then
630 call comm_allreduce(st%st_kpt_mpi_grp, stress_ps_nl)
631 end if
632
633 ! Symmetrize the kinetic stress tensor
634 call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_ps_nl)
635
636 ! Symmetrize the stress tensor if we use k-point symmetries
637 if (hm%kpoints%use_symmetries) then
638 call dsymmetrize_tensor_cart(gr%symm, stress_ps_nl, use_non_symmorphic=.true.)
639 end if
640
641 ! Add the nonlocal energy
642 do idir = 1, ions%space%periodic_dim
643 stress_ps_nl(idir, idir) = stress_ps_nl(idir, idir) + hm%energy%extern_non_local
644 end do
645
646 stress_ps_nl = stress_ps_nl/ions%latt%rcell_volume
647
648 call profiling_out("STRESS_FROM_PSEUDO_NL")
650
651 end subroutine stress_from_pseudo_nonloc
652
653
654 ! -------------------------------------------------------
674 subroutine stress_from_pseudo_local(gr, hm, ions, rho_total, grad_vh, stress_ps_local)
675 type(grid_t), target, intent(in) :: gr
676 type(hamiltonian_elec_t), intent(in) :: hm
677 type(ions_t), intent(in) :: ions
678 real(real64), contiguous, intent(inout) :: rho_total(:)
679 real(real64), intent(in) :: grad_vh(:,:)
680 real(real64), intent(out) :: stress_ps_local(3, 3)
681
682
683 real(real64) :: stress_SR(3, 3), stress_LR(3, 3)
684 real(real64) :: energy_ps_SR, charge, zi
685 real(real64), allocatable :: vloc(:), rvloc(:,:), rho_local_lr(:), rho_lr(:)
686 real(real64), allocatable :: grad_rho(:,:), rho_lr_x(:,:), vlr(:), grad_vlr(:,:)
687 integer :: idir, jdir, iatom
688 type(ps_t), pointer :: spec_ps
689
690 call profiling_in("STRESS_FROM_PSEUDO_LOC")
692
693 ! calculate stress from short-range local pseudopotentials
694 stress_sr = m_zero
695
696 safe_allocate(vloc(1:gr%np))
697 vloc = m_zero
698 safe_allocate(rvloc(1:gr%np, 1:gr%der%dim))
699 rvloc = m_zero
700 do iatom = 1, ions%natoms
701 call epot_local_pseudopotential_sr(gr, ions, iatom, vloc, rvloc)
702 end do
703 safe_deallocate_a(vloc)
704
705 safe_allocate(grad_rho(1:gr%np,1:gr%der%dim))
706 call dderivatives_grad(gr%der, rho_total, grad_rho)
707
708 energy_ps_sr = hm%energy%extern_local
709 do idir = 1, ions%space%periodic_dim
710 do jdir = idir, ions%space%periodic_dim
711 stress_sr(idir, jdir) = stress_sr(idir, jdir) &
712 +dmf_dotp(gr, rvloc(:, jdir), grad_rho(:, idir))
713 end do
714 stress_sr(idir,idir) = stress_sr(idir,idir) + energy_ps_sr
715 end do
716
717 call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_sr)
718
719 stress_sr = stress_sr/ions%latt%rcell_volume
720
721 safe_deallocate_a(rvloc)
722 safe_deallocate_a(grad_rho)
723
724
725 ! calculate stress from long-range local pseudopotentials
726 stress_lr = m_zero
727
728 ! We treat the long-range part of the local potential as the Hartree term
729 ! We first sum the long range densities from atoms
730 safe_allocate(rho_lr(1:gr%np_part))
731 safe_allocate(rho_lr_x(1:gr%np, 1:gr%der%dim))
732 rho_lr = m_zero
733 rho_lr_x = m_zero
734 safe_allocate(rho_local_lr(1:gr%np))
735 do iatom = ions%atoms_dist%start, ions%atoms_dist%end
736 assert(ions%atom(iatom)%species%is_ps())
737 call species_get_long_range_density(ions%atom(iatom)%species, ions%namespace, ions%space, ions%latt, &
738 ions%pos(:, iatom), gr, rho_local_lr, nlr_x=rho_lr_x)
739
740 call lalg_axpy(gr%np, m_one, rho_local_lr, rho_lr)
741 end do
742 safe_deallocate_a(rho_local_lr)
743
744 if (ions%atoms_dist%parallel) then
745 call comm_allreduce(ions%atoms_dist%mpi_grp, rho_lr)
746 call comm_allreduce(ions%atoms_dist%mpi_grp, rho_lr_x)
747 end if
748
749 do idir = 1, ions%space%periodic_dim
750 do jdir = idir, ions%space%periodic_dim
751 stress_lr(idir, jdir) = stress_lr(idir, jdir) + dmf_dotp(gr, rho_lr_x(:,jdir), grad_vh(:, idir))
752 end do
753 end do
754 safe_deallocate_a(rho_lr_x)
755
756 safe_allocate(vlr(1:gr%np_part))
757 call dpoisson_solve(hm%psolver, ions%namespace, vlr, rho_lr, all_nodes = .true.)
758 safe_deallocate_a(rho_lr)
759
760 safe_allocate(grad_vlr(1:gr%np, 1:gr%der%dim))
761 call dderivatives_grad(gr%der, vlr, grad_vlr)
762 safe_deallocate_a(vlr)
763
764 do idir = 1, ions%space%periodic_dim
765 do jdir = idir, ions%space%periodic_dim
766 stress_lr(idir, jdir) = stress_lr(idir, jdir) - dmf_dotp(gr, grad_vh(:,idir), grad_vlr(:, jdir))/m_two/m_pi
767 end do
768 end do
770 call upper_triangular_to_hermitian(ions%space%periodic_dim, stress_lr)
771
772 safe_deallocate_a(grad_vlr)
773
774 ! Contribution from G=0 component of the long-range part
775 !
776 if (ions%space%periodic_dim == 3) then
777 charge = m_zero
778 do iatom = 1, ions%natoms
779 charge = charge + ions%atom(iatom)%species%get_zval()
780 end do
781
782 do iatom = 1, ions%natoms
783 select type(spec => ions%atom(iatom)%species)
784 type is(pseudopotential_t)
785 zi = spec%get_zval()
786 spec_ps => spec%ps
787
788 do idir = 1, ions%space%periodic_dim
789 stress_lr(idir, idir) = stress_lr(idir, idir) &
790 + m_two*m_pi*spec_ps%sigma_erf**2*charge*zi /ions%latt%rcell_volume
791 end do
792 end select
793 end do
794 end if
795
796 stress_lr = stress_lr/ions%latt%rcell_volume
797
798 stress_ps_local = stress_sr + stress_lr
799
800 call profiling_out("STRESS_FROM_PSEUDO_LOC")
802
803 end subroutine stress_from_pseudo_local
804
805 ! -------------------------------------------------------
806 subroutine epot_local_pseudopotential_sr(mesh, ions, iatom, vpsl, rvpsl)
807 class(mesh_t), intent(in) :: mesh
808 type(ions_t), intent(in) :: ions
809 integer, intent(in) :: iatom
810 real(real64), intent(inout) :: vpsl(:)
811 real(real64), intent(inout) :: rvpsl(:,:)
812
813 integer :: ip
814 real(real64) :: radius, vl_ip
815 type(submesh_t) :: sphere
816 type(ps_t), pointer :: ps
817
819
820 if (.not. ions%atom(iatom)%species%is_ps()) then
822 return
823 endif
824
825 call profiling_in("EPOT_LOCAL_PS_SR")
826
827 select type(spec=>ions%atom(iatom)%species)
828 type is(pseudopotential_t)
829
830 ps => spec%ps
831
832 radius = ps%vl%x_threshold*1.05_real64
833
834 call submesh_init(sphere, ions%space, mesh, ions%latt, ions%pos(:, iatom), radius)
835
836 ! Cannot be written (correctly) as a vector expression since for periodic systems,
837 ! there can be values ip, jp such that sphere%map(ip) == sphere%map(jp).
838 do ip = 1, sphere%np
839 vl_ip = spline_eval(ps%vl, sphere%r(ip))
840 vpsl(sphere%map(ip)) = vpsl(sphere%map(ip)) + vl_ip
841 rvpsl(sphere%map(ip), 1:ions%space%periodic_dim) = rvpsl(sphere%map(ip), 1:ions%space%periodic_dim) &
842 + sphere%rel_x(1:ions%space%periodic_dim, ip) * vl_ip
843 end do
844
845 call submesh_end(sphere)
846
847 nullify(ps)
848
849 end select
850
851 call profiling_out("EPOT_LOCAL_PS_SR")
853 end subroutine epot_local_pseudopotential_sr
854
855
856 ! -------------------------------------------------------
871 subroutine stress_from_hubbard(namespace, gr, st, hm, space, rcell_volume, stress_hubbard)
872 type(namespace_t), intent(in) :: namespace
873 type(grid_t), target, intent(in) :: gr
874 type(states_elec_t), intent(inout) :: st
875 type(hamiltonian_elec_t), intent(in) :: hm
876 type(space_t), intent(in) :: space
877 real(real64), intent(in) :: rcell_volume
878 real(real64), intent(out) :: stress_hubbard(3, 3)
879
880 integer :: ik, ist, idir, jdir
881 integer :: ib, minst, maxst
882 type(wfs_elec_t) :: psib, rvu_psib(3), gpsib(3)
883 complex(real64), allocatable :: stress_tmp(:)
884
885 if (hm%lda_u%level == dft_u_none) then
886 stress_hubbard = m_zero
887 return
888 end if
889
890 push_sub_with_profile(stress_from_hubbard)
891
892 assert(st%wfs_type == type_cmplx)
893
894 safe_allocate(stress_tmp(1:st%block_size))
895
896 stress_hubbard = m_zero
897
898 do ik = st%d%kpt%start, st%d%kpt%end
899
900 if (st%kweights(ik) <= m_epsilon) cycle
902 do ib = st%group%block_start, st%group%block_end
903 minst = states_elec_block_min(st, ib)
904 maxst = states_elec_block_max(st, ib)
905
906 call hamiltonian_elec_copy_and_set_phase(hm, gr, st%d%kpt, st%group%psib(ib, ik), psib)
907
908 ! calculate the gradient
909 call zderivatives_batch_grad(gr%der, psib, gpsib, set_bc=.false.)
910
911 ! Get rV_U |\psi> for all atoms
912 do idir = 1, gr%der%dim
913 call psib%copy_to(rvu_psib(idir))
914 call batch_set_zero(rvu_psib(idir))
915 end do
916
917 call zlda_u_rvu(hm%lda_u, gr, space, hm%d, namespace, psib, rvu_psib)
918
919 do idir = 1,3
920 do jdir = idir,3
921 call zmesh_batch_dotp_vector(gr, gpsib(idir), rvu_psib(jdir), stress_tmp)
922
923 do ist = minst, maxst
924 stress_hubbard(idir, jdir) = stress_hubbard(idir, jdir) &
925 + m_two * st%kweights(ik) * st%occ(ist, ik) * real(stress_tmp(ist-minst+1), real64)
926 end do
927
928 end do
929 end do
930
931 do idir = 1, gr%der%dim
932 call rvu_psib(idir)%end()
933 call gpsib(idir)%end()
934 end do
935 call psib%end()
936 end do
937 end do
938
939 safe_deallocate_a(stress_tmp)
940
941 if (st%parallel_in_states .or. st%d%kpt%parallel) then
942 call comm_allreduce(st%st_kpt_mpi_grp, stress_hubbard)
943 end if
944
945 ! Symmetrize the kinetic stress tensor
946 call upper_triangular_to_hermitian(gr%der%dim, stress_hubbard)
947
948 ! Symmetrize the stress tensor if we use k-point symmetries
949 if (hm%kpoints%use_symmetries) then
950 call dsymmetrize_tensor_cart(gr%symm, stress_hubbard)
951 end if
952
953 ! Add the Hubbard energy
954 do idir = 1,3
955 stress_hubbard(idir, idir) = stress_hubbard(idir, idir) + hm%energy%int_dft_u
956 end do
957
958 stress_hubbard = stress_hubbard/rcell_volume
959
960 pop_sub_with_profile(stress_from_hubbard)
961 end subroutine stress_from_hubbard
962
963
964 ! -------------------------------------------------------
965 subroutine output_stress(iunit, space_dim, stress_tensors, all_terms)
966 integer, intent(in) :: iunit
967 integer, intent(in) :: space_dim
968 type(stress_t), intent(in) :: stress_tensors
969 logical, optional, intent(in) :: all_terms
970
971 logical :: write_all_terms
972 character(len=16) :: stress_unit
973
974 write_all_terms = optional_default(all_terms, .true.)
975
976 write(stress_unit, '(4a,i1)') trim(units_abbrev(units_out%energy)), '/', &
977 trim(units_abbrev(units_out%length)), '^', space_dim
978
979 if (mpi_world%is_root()) then
980
981 if (write_all_terms) then
982 write(iunit, '(3a)') 'Kinetic stress tensor [', trim(stress_unit), '] ='
983 call print_stress_tensor(iunit, space_dim, stress_tensors%kinetic)
984 if (space_dim == 3) then
985 write(iunit, '(a, es15.6, 3a)') 'Kinetic pressure sumrule violation: ', &
986 units_from_atomic(units_out%energy, stress_tensors%kinetic_sumrule), &
987 ' [', trim(units_abbrev(units_out%energy)), ']'
988 write(iunit,*)
989 end if
990
991
992 write(iunit, '(3a)') 'Hartree stress tensor [', trim(stress_unit), '] ='
993 call print_stress_tensor(iunit, space_dim, stress_tensors%Hartree)
994 if (space_dim == 3) then
995 write(iunit, '(a, es15.6, 3a)') 'Hartree pressure sumrule violation: ', &
996 units_from_atomic(units_out%energy, stress_tensors%hartree_sumrule), &
997 ' [', trim(units_abbrev(units_out%energy)), ']'
998 write(iunit,*)
999 end if
1000
1001 write(iunit, '(3a)') 'XC stress tensor [', trim(stress_unit), '] ='
1002 call print_stress_tensor(iunit, space_dim, stress_tensors%xc)
1003
1004 write(iunit, '(3a)') 'Local pseudo. stress tensor [', trim(stress_unit), '] ='
1005 call print_stress_tensor(iunit, space_dim, stress_tensors%ps_local)
1006
1007 write(iunit, '(3a)') 'Nonlocal pseudo. stress tensor [', trim(stress_unit), '] ='
1008 call print_stress_tensor(iunit, space_dim, stress_tensors%ps_nl)
1009
1010 write(iunit, '(3a)') 'Ion-ion stress tensor [', trim(stress_unit), '] ='
1011 call print_stress_tensor(iunit, space_dim, stress_tensors%ion_ion)
1012
1013 write(iunit, '(3a)') 'vdW stress tensor [', trim(stress_unit), '] ='
1014 call print_stress_tensor(iunit, space_dim, stress_tensors%vdw)
1015
1016 write(iunit, '(3a)') 'Hubbard stress tensor [', trim(stress_unit), '] ='
1017 call print_stress_tensor(iunit, space_dim, stress_tensors%hubbard)
1018 end if
1019
1020 write(iunit, '(3a)') 'Total stress tensor [', trim(stress_unit), '] ='
1021 call print_stress_tensor(iunit, space_dim, stress_tensors%total)
1022
1023 end if
1024 end subroutine output_stress
1025
1026
1027 subroutine output_pressure(iunit, space_dim, total_stress_tensor)
1028 integer, intent(in) :: iunit
1029 integer, intent(in) :: space_dim
1030 real(real64), intent(in) :: total_stress_tensor(3,3)
1031
1032 integer :: idim
1033 real(real64) :: pressure
1034 character(len=16) :: stress_unit
1035
1036 write(stress_unit, '(4a,i1)') trim(units_abbrev(units_out%energy)), '/', &
1037 trim(units_abbrev(units_out%length)), '^', space_dim
1038
1039 pressure = m_zero
1040 do idim = 1, space_dim
1041 pressure = pressure - total_stress_tensor(idim, idim) / real(space_dim, real64)
1042 end do
1043
1044 write(iunit,'(3a,es16.8)', advance="no") 'Pressure [', trim(stress_unit), '] = ', &
1045 units_from_atomic(units_out%energy/units_out%length**space_dim, pressure)
1046 if (space_dim == 3) then
1047 write(iunit,'(2x,a,f16.8)') 'Pressure [GPa] = ', units_from_atomic(unit_gpa, pressure)
1048 else
1049 write(iunit,*)
1050 end if
1051
1052 end subroutine output_pressure
1053
1054 subroutine print_stress_tensor(ounit, space_dim, tensor)
1055 integer, intent(in) :: ounit
1056 integer, intent(in) :: space_dim
1057 real(real64), intent(in) :: tensor(3,3)
1058
1059 real(real64) :: tensor_with_unit(3,3)
1060 integer :: idim, jdim
1061
1062 tensor_with_unit = units_from_atomic(units_out%energy/units_out%length**space_dim, tensor)
1063
1064 write(ounit,'(a9,2x)', advance="no")"T_{ij}"
1065 do jdim = 1, space_dim
1066 write(ounit,'(i18)', advance="no") jdim
1067 end do
1068 write(ounit,*)
1069 do idim = 1, space_dim
1070 write(ounit,'(i9,2x)', advance="no") idim
1071 do jdim = 1, space_dim
1072 write(ounit,'(es18.9)', advance="no") tensor_with_unit(idim, jdim)
1073 end do
1074 write(ounit,*)
1075 end do
1076 write(ounit,*)
1077
1078 end subroutine print_stress_tensor
1079
1080
1081end module stress_oct_m
1082
1083!! Local Variables:
1084!! mode: f90
1085!! coding: utf-8
1086!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:173
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:188
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:118
subroutine, public batch_set_zero(this, np, async)
fill all mesh functions of the batch with zero
Definition: batch_ops.F90:244
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:124
This module implements a calculator for the density and defines related functions.
Definition: density.F90:122
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
subroutine, public zderivatives_batch_grad(der, ffb, opffb, ghost_update, set_bc, to_cartesian, factor)
apply the gradient to a batch of mesh functions
integer, parameter, public spinors
subroutine, public energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
This subroutine calculates the total energy of the system. Basically, it adds up the KS eigenvalues,...
integer, parameter, public scalar_relativistic_zora
Definition: epot.F90:168
integer, parameter, public fully_relativistic_zora
Definition: epot.F90:168
real(real64), parameter, public m_two
Definition: global.F90:193
real(real64), parameter, public m_zero
Definition: global.F90:191
real(real64), parameter, public m_four
Definition: global.F90:195
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:189
integer, parameter, public independent_particles
Theory level.
Definition: global.F90:237
integer, parameter, public generalized_kohn_sham_dft
Definition: global.F90:237
integer, parameter, public kohn_sham_dft
Definition: global.F90:237
real(real64), parameter, public m_epsilon
Definition: global.F90:207
real(real64), parameter, public m_one
Definition: global.F90:192
This module implements the underlying real-space grid.
Definition: grid.F90:119
subroutine, public hamiltonian_elec_copy_and_set_phase(hm, gr, kpt, psib, psib_with_phase)
Copy a batch to another batch and apply the Bloch phase to it.
This module defines classes and functions for interaction partners.
Definition: io.F90:116
subroutine, public ion_interaction_stress(this, space, latt, atom, natoms, pos, stress_ii)
Computes the contribution to the stress tensor the ion-ion energy.
A module to handle KS potential, without the external potential.
integer, parameter, public dft_u_none
Definition: lda_u.F90:203
subroutine, public zlda_u_rvu(this, mesh, space, d, namespace, psib, gpsib)
This routine computes .
Definition: lda_u.F90:5309
This modules implements the routines for doing constrain DFT for noncollinear magnetism.
integer, parameter, public constrain_none
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:117
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:118
subroutine, public zmesh_batch_dotp_vector(mesh, aa, bb, dot, reduce, cproduct)
calculate the vector of dot-products of mesh functions between two batches
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1091
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
subroutine, public dpoisson_solve(this, namespace, pot, rho, all_nodes, kernel, reset)
Calculates the Poisson equation. Given the density returns the corresponding potential.
Definition: poisson.F90:871
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
Definition: ps.F90:116
subroutine, public species_get_long_range_density(species, namespace, space, latt, pos, mesh, rho, sphere_inout, nlr_x)
subroutine, public species_get_nlcc_grad(species, space, latt, pos, mesh, rho_core_grad, gnlcc_x)
real(real64) function, public spline_eval(spl, x)
Definition: splines.F90:443
This module handles spin dimensions of the states and the k-point distribution.
integer pure function, public states_elec_block_max(st, ib)
return index of last state in block ib
integer pure function, public states_elec_block_min(st, ib)
return index of first state in block ib
This module implements the calculation of the stress tensor.
Definition: stress.F90:120
subroutine stress_from_hartree(gr, space, volume, grad_vh, ehartree, stress_Hartree)
Computes the contribution to the stress tensor from the Hartree energy.
Definition: stress.F90:509
subroutine stress_from_kinetic(gr, space, hm, st, symm, rcell_volume, stress_kin)
Computes the contribution to the stress tensor from the kinetic energy.
Definition: stress.F90:420
subroutine stress_from_pseudo_local(gr, hm, ions, rho_total, grad_vh, stress_ps_local)
Computes the contribution from the local part of the pseudopotential.
Definition: stress.F90:770
subroutine stress_from_xc(energy, rcell_volume, periodic_dim, stress_xc)
Computes the contribution to the stress tensor from the xc energy.
Definition: stress.F90:557
subroutine print_stress_tensor(ounit, space_dim, tensor)
Definition: stress.F90:1150
subroutine, public output_pressure(iunit, space_dim, total_stress_tensor)
Definition: stress.F90:1123
subroutine epot_local_pseudopotential_sr(mesh, ions, iatom, vpsl, rvpsl)
Definition: stress.F90:902
subroutine, public stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
This computes the total stress on the lattice.
Definition: stress.F90:188
subroutine stress_from_hubbard(namespace, gr, st, hm, space, rcell_volume, stress_hubbard)
Computes the contribution to the stress tensor from the Hubbard energy.
Definition: stress.F90:967
subroutine stress_from_xc_nlcc(rcell_volume, gr, st, ions, vxc, stress_xc_nlcc)
Computes the NLCC contribution to the stress tensor from the xc energy.
Definition: stress.F90:587
subroutine stress_from_pseudo_nonloc(gr, st, hm, ions, stress_ps_nl)
Computes the contribution to the stress tensor from the nonlocal part of the pseudopotentials.
Definition: stress.F90:661
subroutine, public output_stress(iunit, space_dim, stress_tensors, all_terms)
Definition: stress.F90:1061
subroutine, public submesh_end(this)
Definition: submesh.F90:738
subroutine, public submesh_init(this, space, mesh, latt, center, rc)
Definition: submesh.F90:282
subroutine, public dsymmetrize_tensor_cart(symm, tensor, use_non_symmorphic)
Symmetric a rank-2 tensor defined in Cartesian space.
type(type_t), public type_cmplx
Definition: types.F90:136
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:225
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
type(unit_t), public unit_gpa
For output pressure in GPa.
Definition: xc.F90:116
logical pure function, public xc_is_energy_functional(xcs)
Is one of the x or c functional is not an energy functional.
Definition: xc.F90:794
pure logical function, public in_family(family, xc_families)
Definition: xc.F90:669
A module that takes care of xc contribution from vdW interactions.
Definition: xc_vdw.F90:118
integer(int64), dimension(5), parameter, public d3_lib_options
VDWCORRECTION options that correspond to the DFT-D3 library.
Definition: xc_vdw.F90:171
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
Describes mesh distribution to nodes.
Definition: mesh.F90:187
A type storing the information and data about a pseudopotential.
Definition: ps.F90:188
The states_elec_t class contains all electronic wave functions.
A submesh is a type of mesh, used for the projectors in the pseudopotentials It contains points on a ...
Definition: submesh.F90:177
batches of electronic states
Definition: wfs_elec.F90:141
int true(void)