From e26be686e8a9be47c292a15422c91a386c787cec Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Thu, 28 Sep 2023 17:12:05 +0200 Subject: [PATCH 1/3] Pressure (#2) * gruneisen | (re-)enable dumping grueneisen parameters on grid * extract_forceconstants | compute compliances, compressibility, bulk modulus - dump compliances to file * compliances | attach to fc and create in %generate * gruneisen_tensors | split out computation of Grueneisen tensors * gruneisen | implement grueneisen tensors, right now without symmetry * gruneisen | implement actual gruneisen tensors, with symmetry * fix | compliances / bulk modulus * phonon_dispersions | store unit cell volume in hdf5 file * extract_forceconstants | fix output for all-electron energies * pressure | add pressure statistics and purely harmonic contribution * pressure | add quasiharmonic/thirdorder term * pressure | symmetrize, outsource write functions * mdsim | use pressure statistics including std. error, remove that from extract_forceconstants * extract_forceconstants | new naming scheme - outfile.new_ucposcar --> outfile.ucposcar.new_refpos for 1. order optimization * solve FC | consecutively fit 1., 2., 3., 4. order - then 1. order is average (symmetrized) force * pressure | compute stress and pressure from force constants - use for relaxation * maint | write 95% error interval (1.96*stderr) * maint | force errors in extract_forceconstants and slight refactoring * pressure | implement (proper) 2nd order (!) and 4th order * maint | remove unused variables, reduce warnings * fix - I don't understand everything * possibly reconcile 38930be and 9072f3e --- src/extract_forceconstants/Makefile | 12 +- src/extract_forceconstants/helper.f90 | 94 +++ src/extract_forceconstants/main.f90 | 614 +++++++++++++++++- src/extract_forceconstants/options.f90 | 6 + src/extract_forceconstants/symmetrize.f90 | 252 +++++++ src/libolle/Makefile | 5 +- src/libolle/ifc_solvers_lsq.f90 | 101 +-- .../type_forceconstant_secondorder.f90 | 12 +- src/libolle/type_forceconstant_thirdorder.f90 | 98 ++- src/libolle/type_mdsim.f90 | 46 +- src/libolle/type_phonon_dispersions.f90 | 53 +- .../type_phonon_dispersions_generation.f90 | 289 ++++++--- src/phonon_dispersion_relations/main.f90 | 4 +- 13 files changed, 1365 insertions(+), 221 deletions(-) create mode 100644 src/extract_forceconstants/helper.f90 create mode 100644 src/extract_forceconstants/symmetrize.f90 diff --git a/src/extract_forceconstants/Makefile b/src/extract_forceconstants/Makefile index 2542bdff..e01f2084 100644 --- a/src/extract_forceconstants/Makefile +++ b/src/extract_forceconstants/Makefile @@ -5,7 +5,9 @@ OBJECT_PATH=../../build/$(CODE)/ OBJS = \ $(OBJECT_PATH)main.o\ -$(OBJECT_PATH)options.o +$(OBJECT_PATH)options.o\ +$(OBJECT_PATH)symmetrize.o\ +$(OBJECT_PATH)helper.o LPATH = -L../../lib $(blaslapackLPATH) $(incLPATHhdf) $(incLPATHmpi) IPATH = -I../../inc/libolle -I../../inc/libflap $(blaslapackIPATH) $(incLPATHhdf) $(incIPATHmpi) @@ -26,7 +28,11 @@ clean: .f90.o: $(F90) $(F90FLAGS) -c $< $(LIBS) -$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o +$(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o $(OBJECT_PATH)symmetrize.o $(OBJECT_PATH)helper.o $(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@ -$(OBJECT_PATH)options.o: +$(OBJECT_PATH)options.o: options.f90 $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ +$(OBJECT_PATH)symmetrize.o: symmetrize.f90 + $(F90) $(F90FLAGS) -c symmetrize.f90 $(LIBS) -o $@ +$(OBJECT_PATH)helper.o: helper.f90 + $(F90) $(F90FLAGS) -c helper.f90 $(LIBS) -o $@ diff --git a/src/extract_forceconstants/helper.f90 b/src/extract_forceconstants/helper.f90 new file mode 100644 index 00000000..bb99411f --- /dev/null +++ b/src/extract_forceconstants/helper.f90 @@ -0,0 +1,94 @@ +!> Tools to symmetrize things, mostly from +!> https://github.com/ollehellman/aims_wrappers/blob/master/src/relax_with_symmetry_constraints/type_irredcell.f90 + +module helper +use konstanter, only: r8, lo_sqtol, lo_exitcode_symmetry, lo_status +use gottochblandat, only: open_file + +implicit none + +contains + +subroutine write_stress_3x3(stress, stress_err, header) + real(r8), dimension(3, 3), intent(in) :: stress, stress_err + character(len=*), intent(in) :: header + character(len=1000) :: opf + + write (*, *) header + opf = "(1X,' average xx xy xz (GPa): ',3(F12.5,' '),' +/- ',3(F12.5,' ') )" + write (*, opf) stress(:, 1), stress_err(:, 1) + opf = "(1X,' average yx yy yz (GPa): ',3(F12.5,' '),' +/- ',3(F12.5,' ') )" + write (*, opf) stress(:, 2), stress_err(:, 2) + opf = "(1X,' average zx zy zz (GPa): ',3(F12.5,' '),' +/- ',3(F12.5,' ') )" + write (*, opf) stress(:, 3), stress_err(:, 3) + write (*, *) +end subroutine + +subroutine write_stress_voigt(stress, stress_std, stress_err, header) + real(r8), dimension(3, 3), intent(in) :: stress, stress_std, stress_err + character(len=*), intent(in) :: header + character(len=1000) :: opf + + write (*, *) header + opf = "(1X,' average xx yy zz xz yz xy (GPa): ',6(F12.5,' '),' +/- ')" + associate (s => stress) + write (*, opf) s(1, 1), s(2, 2), s(3, 3), s(3, 2), s(3, 1), s(2, 1) + end associate + + opf = "(1X,' stddev xx yy zz xz yz xy (GPa): ',6(F12.5,' '))" + associate (s => stress_std) + write (*, opf) s(1, 1), s(2, 2), s(3, 3), s(3, 2), s(3, 1), s(2, 1) + end associate + + opf = "(1X,' stderr xx yy zz xz yz xy (GPa): ',6(F12.5,' '))" + associate (s => stress_err) + write (*, opf) s(1, 1), s(2, 2), s(3, 3), s(3, 2), s(3, 1), s(2, 1) + end associate + + write (*, *) +end subroutine + +subroutine write_file_stress_voigt_time(stress, file, header) + real(r8), dimension(:, :, :), intent(in) :: stress + character(len=*), intent(in) :: header, file + integer :: u, i + + u = open_file('out', file) + write (u, '(a)') header + ! write (u, '(a)') 'timestep,stress_xx,stress_yy,stress_zz,stress_yz,stress_xz,stress_xy' + write (u, '(a)') 'timestep,xx,yy,zz,yz,xz,xy' + associate (s => stress) + do i = 1, size(s, dim=3) + write (u, '(I5,6(",",ES15.8))') i, & + s(1, 1, i), s(2, 2, i), s(3, 3, i), & + s(2, 3, i), s(1, 3, i), s(1, 2, i) + end do + end associate + close (u) + +end subroutine + +subroutine write_file_stress_voigt(stress, stress_std, stress_err, file, header) + real(r8), dimension(3, 3), intent(in) :: stress, stress_std, stress_err + character(len=*), intent(in) :: header, file + integer :: u + + u = open_file('out', file) + write (u, '(a)') header + write (u, '(a)') '# stress (Voigt xx, yy, zz, yz, xz, xy)' + associate (s => stress) + write (u, '(1X,6(1X,ES15.8))') s(1, 1), s(2, 2), s(3, 3), s(3, 2), s(3, 1), s(2, 1) + end associate + write (u, '(a)') '# std. dev. (Voigt xx, yy, zz, yz, xz, xy)' + associate (s => stress_std) + write (u, '(1X,6(1X,ES15.8))') s(1, 1), s(2, 2), s(3, 3), s(3, 2), s(3, 1), s(2, 1) + end associate + write (u, '(a)') '# std. err. (Voigt xx, yy, zz, yz, xz, xy)' + associate (s => stress_err) + write (u, '(1X,6(1X,ES15.8))') s(1, 1), s(2, 2), s(3, 3), s(3, 2), s(3, 1), s(2, 1) + end associate + close (u) + +end subroutine + +end module diff --git a/src/extract_forceconstants/main.f90 b/src/extract_forceconstants/main.f90 index a9406ca3..497da97f 100644 --- a/src/extract_forceconstants/main.f90 +++ b/src/extract_forceconstants/main.f90 @@ -1,10 +1,14 @@ program extract_forceconstants !!{!src/extract_forceconstants/manual.md!}! -use konstanter, only: r8, i8, lo_tol, lo_bohr_to_A, lo_pressure_HartreeBohr_to_GPa, lo_Hartree_to_eV, lo_eV_to_Hartree, & - lo_exitcode_param, lo_force_hartreebohr_to_eVa, lo_iou, lo_forceconstant_1st_HartreeBohr_to_eVA -use gottochblandat, only: tochar, walltime, open_file, lo_chop, lo_does_file_exist, lo_trueNtimes, lo_mean, lo_stddev, & - lo_outerproduct, lo_trace, lo_frobnorm, lo_linear_least_squares +use konstanter, only: r8, i8, lo_iou, lo_tol, lo_bohr_to_A, lo_pressure_HartreeBohr_to_GPa, & + lo_Hartree_to_eV, lo_exitcode_param, & + lo_forceconstant_1st_HartreeBohr_to_eVA, lo_force_hartreebohr_to_eVa, & + lo_kb_Hartree +use gottochblandat, only: tochar, walltime, open_file, lo_chop, lo_does_file_exist, & + lo_trueNtimes, lo_mean, lo_stddev, lo_outerproduct, lo_trace, & + lo_frobnorm, lo_linear_least_squares, lo_invert_real_matrix use mpi_wrappers, only: lo_mpi_helper, lo_stop_gracefully +use lo_memtracker, only: lo_mem_helper use type_crystalstructure, only: lo_crystalstructure use type_forceconstant_firstorder, only: lo_forceconstant_firstorder use type_forceconstant_secondorder, only: lo_forceconstant_secondorder @@ -18,6 +22,9 @@ program extract_forceconstants use lo_memtracker, only: lo_mem_helper use options, only: lo_opts +use symmetrize, only: symmetrize_stress, symmetrize_forces +use helper, only: write_stress_3x3, write_stress_voigt, write_file_stress_voigt_time, & + write_file_stress_voigt use ifc_solvers, only: lo_solve_for_irreducible_ifc, lo_solve_for_irreducible_ifc_fastugly implicit none @@ -312,8 +319,15 @@ program extract_forceconstants ! Now we can get the actual forceconstants from the forcemap and the irreducible representation getfc: block - integer :: i + integer :: i, uu, ii + real(r8) :: compressibility, bulk_modulus real(r8), dimension(3) :: force_norm + real(r8), dimension(6, 6) :: c_ij, s_ij ! elastic constants, compliances + + c_ij = 0.0_r8 + s_ij = 0.0_r8 + compressibility = 0.0_r8 + bulk_modulus = 0.0_r8 if (map%have_fc_singlet) then call map%get_firstorder_forceconstant(uc, fc1) @@ -337,12 +351,38 @@ program extract_forceconstants call map%get_secondorder_forceconstant(uc, fc2, mem, opts%verbosity) if (mw%talk) then call fc2%writetofile(uc, 'outfile.forceconstant') + + ! elastic properties call fc2%get_elastic_constants(uc) + c_ij = fc2%elastic_constants_voigt*lo_pressure_HartreeBohr_to_GPa + s_ij = fc2%elastic_compliances_voigt/lo_pressure_HartreeBohr_to_GPa + write (*, *) '' write (*, *) 'ELASTIC CONSTANTS (GPa):' + do i = 1, 6 - write (*, "(6(3X,F15.5))") lo_chop(fc2%elastic_constants_voigt(:, i)*lo_pressure_HartreeBohr_to_GPa, lo_tol) + write (*, "(6(3X,F15.5))") lo_chop(c_ij(:, i), lo_tol) end do + + ! compressibility and bulk modulus (Nye p. 146) + compressibility = s_ij(1, 1) + s_ij(2, 2) + s_ij(3, 3) + compressibility = compressibility + 2.0_r8*(s_ij(1, 2) + s_ij(2, 3) + s_ij(3, 1)) + bulk_modulus = 1.0_r8/compressibility + write (*, *) '--> Isothermal compressibility (1/GPa):' + write (*, "(3X,F15.5)") compressibility + write (*, *) '--> Isothermal bulk modulus (GPa):' + write (*, "(3X,F15.5)") bulk_modulus + + ! fkdev: write to file + write (*, *) '.. write compliance tensor to outfile.compliance_tensor' + uu = open_file('out', 'outfile.compliance_tensor') + + write (uu, *) "# Compliance tensor in Voigt notation (1/GPa)" + do ii = 1, 6 + write (uu, "(1X,6(1X,F20.15))") s_ij(:, ii) + end do + close (uu) + end if end if if (map%have_fc_triplet) then @@ -623,18 +663,29 @@ program extract_forceconstants update_reference: block ! [ ] worry about long-range later ! [ ] maybe there is a more clever way as usual - integer :: ii, jj + integer :: ii, jj, ip, a1 real(r8), dimension(:), allocatable :: f, s, im + real(r8), dimension(3, sim%na) :: f0, f0_std ! , f2, f3 + real(r8), dimension(3, uc%na) :: f0_prim, f0_prim_sym ! , f2_prim, f2_prim_sym, f3_prim, f3_prim_sym + real(r8), dimension(3, uc%na) :: f0_prim_std, f0_prim_err, f0_prim_err_sym complex(r8), dimension(:, :), allocatable :: D type(lo_qpoint) :: gammapoint + type(lo_crystalstructure) :: uc_new, ss_new if (mw%talk) then write (*, *) '' write (*, *) & - 'Since we have the first order, predict new reference positions.' + '--> since ie have the first order, predict new reference positions', & + ' ind write to outfile.??poscar.new_refpos' end if - associate (n => uc%na) + ! get fresh supercells which we can change + call uc_new%readfromfile('infile.ucposcar') + call ss_new%readfromfile('infile.ssposcar') + call uc_new%classify('wedge', timereversal=.true.) + call ss_new%classify('supercell', uc_new) + + associate (n => uc_new%na) allocate (f(3*n)) allocate (s(3*n)) allocate (im(3*n)) @@ -648,11 +699,8 @@ program extract_forceconstants D = 0.0_r8 gammapoint%r = 0.0_8 - ! make sure supercell knows about unitcell - call ss%classify('supercell', uc) - ! get dynamical matrix - call fc2%dynamicalmatrix(uc, gammapoint, D, mem) + call fc2%dynamicalmatrix(uc_new, gammapoint, D, mem) ! matrix at gamma should be real, sanity check associate (norm => lo_frobnorm(aimag(D))) @@ -664,47 +712,559 @@ program extract_forceconstants end associate ! arrange mass-weighted forces and mass weights - do ii = 1, uc%na - f(3*(ii - 1) + 1:3*ii) = fc1%atom(ii)%m*uc%invsqrtmass(ii) - im(3*(ii - 1) + 1:3*ii) = uc%invsqrtmass(ii) + do ii = 1, uc_new%na + f(3*(ii - 1) + 1:3*ii) = fc1%atom(ii)%m*uc_new%invsqrtmass(ii) + im(3*(ii - 1) + 1:3*ii) = uc_new%invsqrtmass(ii) end do ! predict step as D.T @ s = f call lo_linear_least_squares(transpose(real(D)), f, s) ! undo mass-weighting - s = s*im + s(:) = s*im ! update unitcell positions - do ii = 1, uc%na + do ii = 1, uc_new%na associate (disp => s(3*(ii - 1) + 1:3*ii)) ! in fractional coords - associate (d => uc%cartesian_to_fractional(disp, pbc=.false.), & - r => uc%r(:, ii)) + associate (d => uc_new%cartesian_to_fractional(disp, pbc=.false.), & + r => uc_new%r(:, ii)) r = r - d end associate end associate end do ! update supercell positions - do jj = 1, ss%na - ii = ss%info%index_in_unitcell(jj) + do jj = 1, ss_new%na + ii = ss_new%info%index_in_unitcell(jj) associate (disp => s(3*(ii - 1) + 1:3*ii)) ! add displacement in fractional coords - associate (d => ss%cartesian_to_fractional(disp, pbc=.false.), & - r => ss%r(:, jj)) + associate (d => ss_new%cartesian_to_fractional(disp, pbc=.false.), & + r => ss_new%r(:, jj)) r = r - d end associate end associate - end do + end do ! dump files - call uc%writetofile('outfile.new_ucposcar', 1) - call ss%writetofile('outfile.new_ssposcar', 1) + call uc_new%writetofile('outfile.ucposcar.new_refpos', 1) + call ss_new%writetofile('outfile.ssposcar.new_refpos', 1) + + ! some force statistic + ! I used this to compare to the first-order FC + f0 = 0.0_r8 + f0_std = 0.0_r8 + ! f2 = 0.0_r8 + ! f3 = 0.0_r8 + f0_prim = 0.0_r8 + f0_prim_std = 0.0_r8 + f0_prim_err = 0.0_r8 + f0_prim_err_sym = 0.0_r8 + ! f2_prim = 0.0_r8 + ! f3_prim = 0.0_r8 + do a1 = 1, ss%na + ip = ss%info%index_in_unitcell(a1) + do ii = 1, 3 + f0(ii, a1) = lo_mean(sim%f(ii, a1, :)) + f0_std(ii, a1) = lo_stddev(sim%f(ii, a1, :)) + ! f2(i, a1) = lo_mean(f2t(i, a1, :)) + ! f3(i, a1) = lo_mean(f3t(i, a1, :)) + end do + f0_prim(:, ip) = f0_prim(:, ip) + f0(:, a1) + f0_prim_std(:, ip) = f0_prim(:, ip) + f0(:, a1) + ! f2_prim(:, ip) = f2_prim(:, ip) + f2(:, a1) + ! f3_prim(:, ip) = f3_prim(:, ip) + f3(:, a1) + end do + + f0_prim = f0_prim/ss%na*uc%na + f0_prim_err = f0_prim_std/ss%na*uc%na/sqrt(1.0*sim%nt) + ! f2_prim = f2_prim/ss%na*uc%na + ! f3_prim = f3_prim/ss%na*uc%na + + ! write (lo_iou, *) 'Mean force:' + ! write (*, *) sum(f0_prim, dim=2)/uc%na + ! write (*, *) sum(f2_prim, dim=2)/uc%na + ! write (*, *) sum(f3_prim, dim=2)/uc%na + + ! ! symmetrize the averaged forces + call symmetrize_forces(f0_prim, uc, f0_prim_sym, verbose=.false.) + call symmetrize_forces(f0_prim_err, uc, f0_prim_err_sym, verbose=.false.) + ! call symmetrize_forces(f2_prim, uc, f2_prim_sym) + ! call symmetrize_forces(f3_prim, uc, f3_prim_sym) + + if (mw%talk) then + write (lo_iou, *) + write (lo_iou, *) 'Averaged forces in primitive cell (symmetrized, should be == firstorder FC) with 95% confidence interval:' + do ip = 1, uc%na + write (*, "(1X,I4,3(ES11.3),' +/- ',3(ES11.3))") ip, & + f0_prim_sym(:, ip)*lo_force_hartreebohr_to_eVa, 1.96*f0_prim_err(:, ip)*lo_force_hartreebohr_to_eVa + end do + end if + + ! if (map%have_fc_pair) then + ! write (lo_iou, *) 'Averaged harmonic forces in primitive cell (symmetrized):' + ! do ip = 1, uc%na + ! write (*, *) ip, f2_prim_sym(:, ip)*lo_force_hartreebohr_to_eVa + ! end do + ! end if + + ! if (map%have_fc_triplet) then + ! write (lo_iou, *) 'Averaged third-order forces in primitive cell (symmetrized):' + ! do ip = 1, uc%na + ! write (*, *) ip, f3_prim_sym(:, ip)*lo_force_hartreebohr_to_eVa + ! end do + ! end if end block update_reference end if +! pressure business +if (opts%pressure) then + ! get realspace stress from the force constants: + ! two contributions: + ! - [x] harmonic contribution \Phi_ij u_i u_j + ! - [x] quasi-harmonic contribution \Phi'_ij u_i u_j where \Phi' comes from 3. order FC + + pressure_tdep: block + ! this case is similar to harmonic energy + + type(lo_forceconstant_secondorder) :: fc2_ss + type(lo_forceconstant_thirdorder) :: fc3_ss + type(lo_forceconstant_fourthorder) :: fc4_ss + real(r8), dimension(:, :, :, :), allocatable :: polarfc + real(r8), dimension(3, 3, sim%nt) :: s0t, s2t, spt, s3t, s4t, srt + real(r8), dimension(3, sim%na, sim%nt) :: f0t, fpt, f2t, f3t, f4t + real(r8), dimension(sim%nt) :: e0t, ept, e2t, e3t, e4t, ert, p0t, p2t, ppt, p3t, p4t, prt + real(r8), dimension(3, 3, 3) :: m3 + real(r8), dimension(3, 3, 3, 3) :: m4 + real(r8), dimension(3, 3) :: m2, s0, s2, sp, s3, s4, sr + real(r8), dimension(3, 3) :: s0_std, s2_std, s3_std, s4_std, sr_std + real(r8), dimension(3, 3) :: s0_err, s2_err, s3_err, s4_err, sr_err + real(r8), dimension(3, 3) :: s0_sym, s2_sym, s3_sym, sr_sym ! s4_sym + real(r8), dimension(3, 3) :: s0_std_sym, s2_std_sym, s3_std_sym, sr_std_sym ! s4_std_sym + real(r8), dimension(3, 3) :: s0_err_sym, s2_err_sym, s3_err_sym, sr_err_sym ! s4_err_sym + real(r8), dimension(3) :: v0, u1, u2, u3, u4, rv1, rv2, rv3, rv4 + !> averaged force + ! real(r8), dimension(3, sim%na) :: f0, f2, f3, f4, f0_std + ! real(r8), dimension(3, uc%na) :: f0_prim, f0_prim_sym, f2_prim, f2_prim_sym, f3_prim, f3_prim_sym + ! real(r8), dimension(3, uc%na) :: f0_prim_std, f0_prim_err + ! real(r8) :: dfp, df2, df3, df4 + real(r8) :: n_samples, energy, baseline, tomev, toev, volume + real(r8) :: p0, p2, pp, p3, p4, pr + real(r8) :: p0_std, p2_std, pp_std, p3_std, p4_std, pr_std + real(r8) :: p0_err, p2_err, pp_err, p3_err, p4_err, pr_err + integer :: i, j, t, a1, a2, a3, a4, ii, i1, i2, i3, i4, i5, u, ctr, ctrtot + + ! Make sure we have the supercell and MD data + if (ss%na .lt. 0) then + call ss%readfromfile('infile.ssposcar') + call ss%classify('supercell', uc) + end if + if (sim%na .lt. 0) then + call sim%read_from_file(verbosity=2, stride=opts%stride) + call sim%remove_force_and_center_of_mass_drift() + end if + + ! system volume and number of samples + n_samples = 1.0_r8*sim%nt + volume = ss%volume + + ! Remap the forceconstants to the supercell + if (map%have_fc_pair) call fc2%remap(uc, ss, fc2_ss) + ! fkdev: thirdorder will become sth like secondorder, later problem + if (map%have_fc_triplet) call fc3%remap(uc, ss, fc3_ss) + if (map%have_fc_quartet) call fc4%remap(uc, ss, fc4_ss) + if (map%polar .gt. 0) then + allocate (polarfc(3, 3, ss%na, ss%na)) + polarfc = 0.0_r8 + call fc2%supercell_longrange_dynamical_matrix_at_gamma(ss, polarfc, 1E-15_r8) + end if + + ! init + p0t = 0.; p2t = 0.; ppt = 0.; p3t = 0.; p4t = 0.; prt = 0.; + p0_std = 0.; p2_std = 0.; pp_std = 0.; p3_std = 0.; p0_err = 0.; p2_err = 0.; pp_err = 0.; p3_err = 0.; + s0t = 0.; spt = 0.; s2t = 0.; s3t = 0.; s4t = 0.; sr = 0.; srt = 0.; + s0_std = 0.; s2_std = 0.; s3_std = 0.; sr_std = 0.; + s0_err = 0.; s2_err = 0.; s3_err = 0.; sr_err = 0.; + s0_sym = 0.; s2_sym = 0.; s3_sym = 0.; sr_sym = 0.; + s0_std_sym = 0.; s2_std_sym = 0.; s3_std_sym = 0.; + s0_err_sym = 0.; s2_err_sym = 0.; s3_err_sym = 0.; + f0t = 0.; fpt = 0.; f2t = 0.; f3t = 0.; f4t = 0.; + e0t = 0.; ept = 0.; e2t = 0.; e3t = 0.; e4t = 0.; ert = 0.; + p0 = 0.; p2 = 0.; pp = 0.; p3 = 0.; p4 = 0.; p0 = 0.; + ! Unit conversion + tomev = lo_Hartree_to_eV*1000/ss%na + toev = lo_Hartree_to_eV/ss%na + + ctr = 0 + ctrtot = 0 + do t = 1, sim%nt + if (mod(t, mw%n) .ne. mw%r) cycle + ctrtot = ctrtot + 1 + end do + + ! Calculate forces and stresses + do t = 1, sim%nt + ! make it parallel to not confuse anyone + if (mod(t, mw%n) .ne. mw%r) cycle + ! Copy of DFT energy, force, stress + e0t(t) = sim%stat%potential_energy(t) + f0t(:, :, t) = sim%f(:, :, t) + s0t(:, :, t) = sim%stat%stress(:, :, t) + ! then the harmonic terms + if (map%have_fc_pair) then + energy = 0.0_r8 + s2 = 0.0_r8 + do a1 = 1, ss%na + v0 = 0.0_r8 + u1 = sim%u(:, a1, t) + do ii = 1, fc2_ss%atom(a1)%n + a2 = fc2_ss%atom(a1)%pair(ii)%i2 + m2 = fc2_ss%atom(a1)%pair(ii)%m + rv2 = fc2_ss%atom(a1)%pair(ii)%r + v0 = v0 - matmul(m2, sim%u(:, a2, t)) + ! fkdev: used fixed reference position + do i1 = 1, 3 + do i2 = 1, 3 + do i3 = 1, 3 + associate ( & + s => s2(i1, i2), & + prefactor => 1.0_r8/volume, & + phi_prime => m2(i3, i1)*u1(i3)*rv2(i2) & + ) + s = s + prefactor*phi_prime + end associate + end do + end do + end do + + end do + associate (f => v0) + energy = energy - 0.5_r8*dot_product(u1, f) + f2t(:, a1, t) = f + s2 = s2 - lo_outerproduct(u1, f)/volume + end associate + end do + e2t(t) = energy + s2t(:, :, t) = s2 + end if + ! Possible polar term? + if (fc2%polar) then + energy = 0.0_r8 + sp = 0.0_r8 + do a1 = 1, ss%na + v0 = 0.0_r8 + do a2 = 1, ss%na + v0 = v0 - matmul(polarfc(1:3, 1:3, a1, a2), sim%u(1:3, a2, t)) + end do + associate (u => sim%u(:, a1, t), f => v0) + energy = energy - dot_product(u, v0)*0.5_r8 + fpt(:, a1, t) = f + sp = sp - 0.5_r8/volume*(lo_outerproduct(u, f) + lo_outerproduct(f, u)) + end associate + + end do + ept(t) = energy + spt(:, :, t) = sp + end if + ! triplet term gives quasi-harmonic contribution + if (map%have_fc_triplet) then + energy = 0.0_r8 + s3 = 0.0_r8 + do a1 = 1, fc3_ss%na + v0 = 0.0_r8 + do i = 1, fc3_ss%atom(a1)%n + m3 = fc3_ss%atom(a1)%triplet(i)%m + a2 = fc3_ss%atom(a1)%triplet(i)%i2 + a3 = fc3_ss%atom(a1)%triplet(i)%i3 + u1 = sim%u(:, a1, t) + u2 = sim%u(:, a2, t) + u3 = sim%u(:, a3, t) + + ! positions + rv1 = fc3_ss%atom(a1)%triplet(i)%rv1 + rv2 = fc3_ss%atom(a1)%triplet(i)%rv2 + rv3 = fc3_ss%atom(a1)%triplet(i)%rv3 + + do ii = 1, 3 + do i2 = 1, 3 + do i3 = 1, 3 + v0(ii) = v0(ii) - m3(ii, i2, i3)*u2(i2)*u3(i3) + do i4 = 1, 3 + associate ( & + s => s3(ii, i2), & + prefactor => 0.5_r8/volume, & + ! fkdev: use r_i = r0_i + u_i instead of r0_i + phi_prime1 => m3(i3, i4, i2)*u1(i3)*u2(i4)*rv3(ii) & + ! phi_prime2 => m3(i3, i2, i4)*u1(i3)*u3(i4)*rv2(i1) & + ) + ! s = s + prefactor*(phi_prime1 + phi_prime2 )/2.0_r8 + s = s + prefactor*phi_prime1 + end associate + end do + end do + end do + end do + end do + v0 = v0*0.5_r8 + f3t(:, a1, t) = v0 + + ! fkdev: add energy term to stress similar to 2nd order + ! comment: doesn't seem to reduce variance/fluctuations + associate (u => sim%u(:, a1, t), f => v0) + s3 = s3 - lo_outerproduct(u, f)/volume + end associate + energy = energy - dot_product(v0, sim%u(:, a1, t))/3.0_r8 + ! fkdev: I think this cannot work b.c. p.b.c but I tried + ! sft(:, :, t) = sft(:, :, t) - lo_outerproduct(ss%rcart(:, a1), v0)/volume/2.0_r8 + end do + e3t(t) = energy + s3t(:, :, t) = s3 + end if + ! quartet term + if (map%have_fc_quartet) then + energy = 0.0_r8 + s4 = 0.0_r8 + do a1 = 1, fc4_ss%na + v0 = 0.0_r8 + do i = 1, fc4_ss%atom(a1)%n + m4 = fc4_ss%atom(a1)%quartet(i)%m + a2 = fc4_ss%atom(a1)%quartet(i)%i2 + a3 = fc4_ss%atom(a1)%quartet(i)%i3 + a4 = fc4_ss%atom(a1)%quartet(i)%i4 + u1 = sim%u(:, a1, t) + u2 = sim%u(:, a2, t) + u3 = sim%u(:, a3, t) + u4 = sim%u(:, a4, t) + + ! positions + rv3 = fc4_ss%atom(a1)%quartet(i)%rv3 + rv4 = fc4_ss%atom(a1)%quartet(i)%rv4 + + do ii = 1, 3 + do i2 = 1, 3 + do i3 = 1, 3 + do i4 = 1, 3 + v0(ii) = v0(ii) - m4(ii, i2, i3, i4)*u2(i2)*u3(i3)*u4(i4) + do i5 = 1, 3 + associate ( & + s => s4(ii, i2), & + prefactor => 1.0_r8/volume/6.0_r8, & + phi_prime => u1(i3)*u2(i4)*u3(i5)*m4(i3, i4, i5, ii)*rv4(i2) & + ) + s = s + prefactor*phi_prime + end associate + + end do + end do + end do + end do + end do + end do + v0 = v0/6.0_r8 + f4t(:, a1, t) = v0 + + ! fkdev: add energy term to stress similar to 2nd order + ! associate (u => sim%u(:, a1, t), f => v0) + ! s4 = s4 - lo_outerproduct(u, f)/volume + ! end associate + + energy = energy - dot_product(v0, sim%u(:, a1, t))/4.0_r8 + end do + e4t(t) = energy + s4t(:, :, t) = s4 + end if + end do + + ! sync across ranks + call mw%allreduce('sum', f0t) + call mw%allreduce('sum', fpt) + call mw%allreduce('sum', f2t) + call mw%allreduce('sum', f3t) + call mw%allreduce('sum', f4t) + call mw%allreduce('sum', e0t) + call mw%allreduce('sum', ept) + call mw%allreduce('sum', e2t) + call mw%allreduce('sum', e3t) + call mw%allreduce('sum', e4t) + call mw%allreduce('sum', s0t) + call mw%allreduce('sum', s2t) + call mw%allreduce('sum', spt) + call mw%allreduce('sum', s3t) + call mw%allreduce('sum', s4t) + + ! Subtract a baseline to get sensible numbers to work with + ert = e0t - e2t - e3t - e4t - ept + ! stress residual + srt = s0t - s2t - spt - s3t - s4t + baseline = lo_mean(ert) + e0t = e0t - baseline + + ! do everything else in GPa for digestive purposes + s0t = s0t*lo_pressure_HartreeBohr_to_GPa + s2t = s2t*lo_pressure_HartreeBohr_to_GPa + spt = spt*lo_pressure_HartreeBohr_to_GPa + s3t = s3t*lo_pressure_HartreeBohr_to_GPa + s4t = s4t*lo_pressure_HartreeBohr_to_GPa + srt = srt*lo_pressure_HartreeBohr_to_GPa + + do i = 1, 3 + do j = 1, 3 + s0(i, j) = lo_mean(s0t(i, j, :)) + s2(i, j) = lo_mean(s2t(i, j, :)) + s3(i, j) = lo_mean(s3t(i, j, :)) + s4(i, j) = lo_mean(s4t(i, j, :)) + sr(i, j) = lo_mean(srt(i, j, :)) + s0_std(i, j) = lo_stddev(s0t(i, j, :)) ! /sqrt(n_samples) + s2_std(i, j) = lo_stddev(s2t(i, j, :)) ! /sqrt(n_samples) + s3_std(i, j) = lo_stddev(s3t(i, j, :)) ! /sqrt(n_samples) + s4_std(i, j) = lo_stddev(s4t(i, j, :)) + sr_std(i, j) = lo_stddev(srt(i, j, :)) ! /sqrt(n_samples) + end do + p0t = p0t - s0t(i, i, :)/3.0_r8 + p2t = p2t - s2t(i, i, :)/3.0_r8 + ppt = ppt - spt(i, i, :)/3.0_r8 + p3t = p3t - s3t(i, i, :)/3.0_r8 + p4t = p4t - s4t(i, i, :)/3.0_r8 + prt = prt - srt(i, i, :)/3.0_r8 + end do + + ! error + s0_err = s0_std/sqrt(n_samples) + s2_err = s2_std/sqrt(n_samples) + s3_err = s3_std/sqrt(n_samples) + s4_err = s4_std/sqrt(n_samples) + sr_err = sr_std/sqrt(n_samples) + + ! pressure + p0 = lo_mean(p0t) + p2 = lo_mean(p2t) + pp = lo_mean(ppt) + p3 = lo_mean(p3t) + p4 = lo_mean(p4t) + pr = lo_mean(prt) + p0_std = lo_stddev(p0t) + p2_std = lo_stddev(p2t) + pp_std = lo_stddev(ppt) + p3_std = lo_stddev(p3t) + p4_std = lo_stddev(p4t) + pr_std = lo_stddev(prt) + p0_err = p0_std/sqrt(n_samples) + p2_err = p2_std/sqrt(n_samples) + pp_err = pp_std/sqrt(n_samples) + p3_err = p3_std/sqrt(n_samples) + p4_err = p4_std/sqrt(n_samples) + pr_err = pr_std/sqrt(n_samples) + + ! symmetrize things + call symmetrize_stress(s0, uc, s0_sym) + call symmetrize_stress(s2, uc, s2_sym) + call symmetrize_stress(s3, uc, s3_sym) + call symmetrize_stress(sr, uc, sr_sym) + call symmetrize_stress(s0_std, uc, s0_std_sym) + call symmetrize_stress(s2_std, uc, s2_std_sym) + call symmetrize_stress(s3_std, uc, s3_std_sym) + call symmetrize_stress(sr_std, uc, sr_std_sym) + call symmetrize_stress(s0_err, uc, s0_err_sym) + call symmetrize_stress(s2_err, uc, s2_err_sym) + call symmetrize_stress(s3_err, uc, s3_err_sym) + call symmetrize_stress(sr_err, uc, sr_err_sym) + + ! Now I can report a little? + if (mw%talk) then + write (*, *) '' + write (*, "(1X,A,11X,A,10X,A,7X,A)") 'ENERGIES:', 'rms', 'stddev', 'stddev(residual) (meV/atom)' + write (*, '(1X,A15,2(1X,F12.6),5X,A)') ' input:', sqrt(lo_mean(e0t**2))*tomev, lo_stddev(e0t)*tomev, '-' + write (*, '(1X,A15,3(1X,F12.6))') 'second order:', sqrt(lo_mean(e2t**2))*tomev, lo_stddev(e2t)*tomev, lo_stddev(e0t - e2t)*tomev + write (*, '(1X,A15,3(1X,F12.6))') ' polar:', sqrt(lo_mean(ept**2))*tomev, lo_stddev(ept)*tomev, lo_stddev(e0t - ept - e2t)*tomev + write (*, '(1X,A15,3(1X,F12.6))') 'third order:', sqrt(lo_mean(e3t**2))*tomev, lo_stddev(e3t)*tomev, lo_stddev(e0t - ept - e2t - e3t)*tomev + write (*, '(1X,A15,3(1X,F12.6))') 'fourth order:', sqrt(lo_mean(e4t**2))*tomev, lo_stddev(e4t)*tomev, lo_stddev(e0t - ept - e2t - e3t - e4t)*tomev + + write (*, "(1X,A,11X,A,10X,A,5X,A)") 'PRESSURE: ', 'mean', 'stddev', 'stdev (residual) (GPa)' + write (*, '(1X,A15,2(1X,F12.6),5X,A)') ' input:', p0, p0_std, '-' + write (*, '(1X,A15,3(1X,F12.6),5X,A)') '2nd order:', p2, p2_std, lo_stddev(p0t - p2t) + write (*, '(1X,A15,3(1X,F12.6),5X,A)') ' polar:', pp, pp_std, lo_stddev(p0t - p2t - ppt) + write (*, '(1X,A15,3(1X,F12.6),5X,A)') '3rd order:', p3, p3_std, lo_stddev(p0t - p2t - ppt - p3t) + write (*, '(1X,A15,3(1X,F12.6),5X,A)') '4th order:', p4, p4_std, lo_stddev(p0t - p2t - ppt - p3t - p4t) + write (*, '(1X,A15,2(1X,F12.6),5X,A)') ' residual:', pr, pr_std, '-' + + write (*, *) + associate (prefactor => 2.0_r8/sim%na/lo_kb_Hartree/3.0_r8) + write (*, "(1X,A,F12.6,A,F12.6,A)") 'Harmonic temperature: ', & + prefactor*lo_mean(e2t), ' +/- ', prefactor*lo_stddev(e2t)/sqrt(n_samples), ' K' + end associate + + write (*, *) + associate (prefactor => 2.0_r8/3.0_r8/volume*lo_pressure_HartreeBohr_to_GPa) + write (*, "(1X,A)") 'Harmonic pressure from energy p = 2/3 E_pot / V: ' + write (*, "(1X,F12.6,A,F12.6,A)") & + prefactor*lo_mean(e2t), ' +/- ', prefactor*lo_stddev(e2t)/sqrt(n_samples), ' GPa' + end associate + write (*, *) + + ! call write_stress_3x3(s0, ds0, 'DFT stress') + ! call write_stress_3x3(s0_sym, ds0_sym, 'DFT stress (symmetrized') + call write_stress_voigt(s0, s0_std, s0_err, 'DFT stress (GPa)') + call write_stress_voigt(s0_sym, s0_std_sym, s0_err_sym, 'DFT stress (symmetrized)') + ! call write_stress_voigt(s2, s2_std, s2_err, 'Harmonic FC stress (2nd order)') + call write_stress_voigt(s2_sym, s2_std_sym, s2_err_sym, 'Harmonic FC stress (2nd order) (symmetrized)') + if (map%have_fc_triplet) then + ! call write_stress_voigt(s3, s3_std, s3_err, 'Anharmonic FC stress (3rd order)') + call write_stress_voigt(s3_sym, s3_std_sym, s3_err_sym, 'Anharmonic FC stress (3rd order) (symetrized)') + ! call write_stress_voigt(sr, sr_std, sr_err, 'Residual stress') + call write_stress_voigt(sr_sym, sr_std_sym, sr_err_sym, 'Residual stress (symmetrized)') + end if + + ! dump pressure to file + write (lo_iou, *) '... dump pressure to file' + u = open_file('out', 'outfile.pressure.csv') + write (u, '(a)') '# unit: GPa' + write (u, '(a)') 'timestep,pressure_dft,pressure_fc2,pressure_fc3,pressure_res,pressure_fc4' + do i = 1, sim%nt + write (u, '(I5,5(",",E15.8))') i, p0t(i), p2t(i), p3t(i), prt(i), p4t(i) + end do + close (u) + + ! dump stress to file for each time step + write (lo_iou, *) '... dump DFT stress to file' + call write_file_stress_voigt_time(s0t, 'outfile.stress_dft.csv', '# DFT stress, unit: GPa') + write (lo_iou, *) '... dump FC2 stress to file' + call write_file_stress_voigt_time(s2t, 'outfile.stress_fc2.csv', '# FC2 stress, unit: GPa') + if (map%have_fc_triplet) then + ! call write_stress_voigt(s3, s3_std, s3_err, 'Anharmonic FC stress (3rd order)') + write (lo_iou, *) '... dump FC3 stress to file' + call write_file_stress_voigt_time(s3t, 'outfile.stress_fc3.csv', '# FC3 stress, unit: GPa') + write (lo_iou, *) '... dump res. stress to file' + call write_file_stress_voigt_time(srt, 'outfile.stress_res.csv', '# Residual stress (DFT-QHA), unit: GPa') + end if + if (map%have_fc_quartet) then + write (lo_iou, *) '... dump FC4 stress to file' + call write_file_stress_voigt_time(s4t, 'outfile.stress_fc4.csv', '# FC4 stress, unit: GPa') + + end if + + ! dump symmetrized stress for relaxation incl. std. dev. and std. err. + ! DFT stress + write (lo_iou, *) '... dump symmetrized DFT stress to file' + call write_file_stress_voigt(s0_sym, s0_std_sym, s0_err_sym, 'outfile.stress_dft', '# DFT stress (symmetrized), unit: GPa') + ! 2nd order FC stress + write (lo_iou, *) '... dump symmetrized FC2 stress to file' + call write_file_stress_voigt(s2_sym, s2_std_sym, s2_err_sym, 'outfile.stress_fc2', '# FC2 stress (symmetrized), unit: GPa') + if (map%have_fc_triplet) then + ! 3rd order FC stress + write (lo_iou, *) '... dump symmetrized FC3 stress to file' + call write_file_stress_voigt(s3_sym, s3_std_sym, s3_err_sym, 'outfile.stress_fc3', '# FC3 stress (symmetrized), unit: GPa') + ! residual stress (DFT - FC2 - FC3) + write (lo_iou, *) '... dump symmetrized residual stress to file' + call write_file_stress_voigt(sr_sym, sr_std_sym, sr_err_sym, 'outfile.stress_res', '# Residual stress (DFT-FC) (symmetrized), unit: GPa') + end if + + end if + + end block pressure_tdep + +end if + if (mw%talk) then write (*, *) '' write (*, *) 'Done in ', tochar(walltime() - t0), 's' diff --git a/src/extract_forceconstants/options.f90 b/src/extract_forceconstants/options.f90 index 3d89a77e..6f75241a 100644 --- a/src/extract_forceconstants/options.f90 +++ b/src/extract_forceconstants/options.f90 @@ -68,6 +68,8 @@ module options logical :: fake_dielectric = .false. ! developer mode logical :: devmode = .false. + ! extract realspace pressure from force constants + logical :: pressure = .false. contains procedure :: parse @@ -247,6 +249,9 @@ subroutine parse(opts) if (lo_status .ne. 0) stop call cli%add(switch='--developermode', switch_ab='-dev', hidden=.true., help='dev. mode', & required=.false., act='store_true', def='.false.', error=lo_status) + call cli%add(switch='--pressure', & + help='Extract TDEP pressure from samples for lattice expansion etc.', & + required=.false., act='store_true', def='.false.', error=lo_status) if (lo_status .ne. 0) stop ! actually parse it @@ -306,6 +311,7 @@ subroutine parse(opts) call cli%get(switch='--temperature', val=opts%temperature) call cli%get(switch='--fakediel', val=opts%fake_dielectric) call cli%get(switch='--developermode', val=opts%devmode) + call cli%get(switch='--pressure', val=opts%pressure) if (lo_status .ne. 0) stop diff --git a/src/extract_forceconstants/symmetrize.f90 b/src/extract_forceconstants/symmetrize.f90 new file mode 100644 index 00000000..484f3c84 --- /dev/null +++ b/src/extract_forceconstants/symmetrize.f90 @@ -0,0 +1,252 @@ +!> Tools to symmetrize things, mostly from +!> https://github.com/ollehellman/aims_wrappers/blob/master/src/relax_with_symmetry_constraints/type_irredcell.f90 + +module symmetrize + +use konstanter, only: r8, lo_sqtol, lo_exitcode_symmetry, lo_status +use type_crystalstructure, only: lo_crystalstructure +use gottochblandat, only: lo_identitymatrix, lo_nullspace_coefficient_matrix, & + lo_stop_gracefully, lo_unflatten_2tensor, lo_chop, tochar, & + lo_flattentensor +use type_blas_lapack_wrappers, only: lo_dgels +implicit none + +contains + +!> symmetrize 3x3 stress/strain +subroutine symmetrize_stress(stress, uc, stress_symmetric, verbose) + !> the matrix that should be symmetrized + real(r8), dimension(3, 3), intent(in) :: stress + real(r8), dimension(3, 3), intent(out) :: stress_symmetric + logical, intent(in), optional :: verbose + logical :: talk + !> the crystal structure with symmetry + type(lo_crystalstructure), intent(in) :: uc + real(r8), dimension(:, :), allocatable :: coeffM, wM + real(r8), dimension(9, 1) :: wV + + real(r8), dimension(:), allocatable :: stress_irred_solution + real(r8), dimension(9, 9) :: invarM, I9, T, rotM + real(r8), dimension(3, 3) :: m0 + real(r8), dimension(9) :: dv + integer :: i, j, n_irred + + talk = .false. + if (present(verbose)) then + if (verbose) talk = .true. + end if + + ! Transposition operator + T(:, 1) = [1, 0, 0, 0, 0, 0, 0, 0, 0] + T(:, 2) = [0, 0, 0, 1, 0, 0, 0, 0, 0] + T(:, 3) = [0, 0, 0, 0, 0, 0, 1, 0, 0] + T(:, 4) = [0, 1, 0, 0, 0, 0, 0, 0, 0] + T(:, 5) = [0, 0, 0, 0, 1, 0, 0, 0, 0] + T(:, 6) = [0, 0, 0, 0, 0, 0, 0, 1, 0] + T(:, 7) = [0, 0, 1, 0, 0, 0, 0, 0, 0] + T(:, 8) = [0, 0, 0, 0, 0, 1, 0, 0, 0] + T(:, 9) = [0, 0, 0, 0, 0, 0, 0, 0, 1] + ! identity matrix + call lo_identitymatrix(I9) + ! now add all the operations of the lattice + invarM = 0.0_r8 + do i = 1, uc%sym%n + rotM = expandoperation_pair(uc%sym%op(i)%m) + invarM = invarM + rotM - I9 + end do + ! and transposition symmetry + invarM = invarM + T - I9 + ! and project out the nullspace + call lo_nullspace_coefficient_matrix(invarM, coeffM, n_irred, tolerance=lo_sqtol) + + ! It's weird if it's not at least 1, or at most 6, since the number of zero singular values + ! should be the number of degrees of freedom. + if (n_irred .lt. 1 .or. n_irred .gt. 6) then + call lo_stop_gracefully(['Strain tensor must have at least 1 and at most 6 degrees of freedom'], & + lo_exitcode_symmetry, __FILE__, __LINE__) + end if + ! and a little space for the solution + allocate (stress_irred_solution(n_irred)) + stress_irred_solution = 0.0_r8 + + if (talk) then + write (*, *) '... found ', tochar(n_irred), ' degrees of freedom for the lattice:' + do i = 1, n_irred + write (*, *) 'x'//tochar(i)//':' + dv = 0.0_r8 + dv(i) = 1.0_r8 + m0 = lo_unflatten_2tensor(matmul(coeffM, dv(1:n_irred))) + do j = 1, 3 + write (*, "(3(2X,F16.9))") m0(:, j) + end do + end do + end if + + ! now solve for the irreducible components + ! Solve for the irreducible components + allocate (wM(9, n_irred)) + wM(:, :) = coeffM + wV(:, 1) = lo_flattentensor(stress) + call lo_dgels(wM, wV, info=lo_status) + if (lo_status .ne. 0) then + write (*, *) 'Lapack failed for the lattice' + stop + end if + stress_irred_solution(:) = wV(1:n_irred, 1) + ! Store the symetrized stress + stress_symmetric(1:3, 1:3) = lo_chop(lo_unflatten_2tensor(matmul(coeffM, stress_irred_solution)), 1E-12_r8) +end subroutine + +!> symmetrize 3xN forces in unitcell +subroutine symmetrize_forces(forces, uc, forces_sym, verbose) + !> the matrix that should be symmetrized + real(r8), dimension(:, :), intent(in) :: forces + real(r8), dimension(:, :), intent(out) :: forces_sym + logical, intent(in), optional :: verbose + logical :: talk + !> the crystal structure with symmetry + type(lo_crystalstructure), intent(in) :: uc + real(r8), dimension(:, :), allocatable :: coeffM, dr + + real(r8), dimension(:), allocatable :: forces_irred + real(r8), dimension(3*uc%na, 3*uc%na) :: invarM, Id, rotM + real(r8), dimension(:), allocatable :: dv + integer :: i, j, o, n_irred, dof + + if ((present(verbose)) .and. (verbose)) then + talk = .true. + else + talk = .false. + end if + + call lo_identitymatrix(Id) + + ! Build a matrix representation how all the internal positions transform under an + ! operation as large block matrices: + invarM = 0.0_r8 + do o = 1, uc%sym%n + rotM = 0.0_r8 + do i = 1, uc%na + ! the semimagic %fmap here keeps track of which atom is transform to which atom + ! under each operation. I can mimic this switcharoo by putting the small 3x3 + ! rotation matrix in the right place in the big matrix. That way I don't have to + ! deal with translations explicitly, and life becomes much easier. + j = uc%sym%op(o)%fmap(i) + !rotM((j-1)*3+1:j*3,(i-1)*3+1:3*i)=transpose(p%sym%op(o)%fm) + ! Mathematica confused me again, this is how it's supposed to be. Hmm. + rotM((i - 1)*3 + 1:i*3, (j - 1)*3 + 1:3*j) = transpose(uc%sym%op(o)%fm) + end do + invarM = invarM + rotM - Id + end do + + ! Then add the zero-sum constraint? You should convince yourself that this really is the + ! matrix representation of adding all forces together. + rotM = 0.0_r8 + do i = 1, uc%na + do j = 1, 3 + rotM((i - 1)*3 + j, j) = 1.08 + end do + end do + invarM = invarM + rotM + + ! rotational invariance? nobody knows + + ! Now we can project out the nullspace of this matrix + call lo_nullspace_coefficient_matrix(invarM, coeffM, n_irred, tolerance=lo_sqtol) + + if (talk) then + write (*, *) '... found ', tochar(n_irred), ' internal degrees of freedom' + dof = uc%na*3 + allocate (dv(dof)) + allocate (dr(3, uc%na)) + do i = 1, n_irred + write (*, *) 'x'//tochar(i)//': (fractional, Cartesian)' + dv = 0.0_r8 + dv(i) = 1.0_r8 + dv = matmul(coeffM, dv(1:n_irred)) + do j = 1, uc%na + dr(:, j) = dv((j - 1)*3 + 1:j*3) + write (*, "(6(2X,F15.9))") dr(:, j), uc%fractional_to_cartesian(dr(:, j)) + end do + end do + end if + + ! now solve for the irreducible components if there is anything to solve + if (n_irred .gt. 0) then + fixforces: block + real(r8), dimension(:, :), allocatable :: mA, mB + real(r8), dimension(:), allocatable :: vB + real(r8), dimension(3) :: v0 + integer :: l, n + + n = size(coeffM, 1) + if (size(forces, 1)*size(forces, 2) .ne. n) then + write (*, *) 'Wrong number of forces' + stop + end if + + ! and some space for a solution + allocate (forces_irred(n_irred)) + forces_irred = 0.0_r8 + + ! normal least squares for forces in fractional coords + allocate (mA(n, n_irred)) + allocate (mB(n, 1)) + l = 0 + do i = 1, size(forces, 2) + associate (f_frac => matmul(uc%inv_latticevectors, forces(:, i))) + do j = 1, 3 + l = l + 1 + mB(l, 1) = f_frac(j) + end do + end associate + end do + mA(:, :) = coeffM + call lo_dgels(mA, mB, info=lo_status) + if (lo_status .ne. 0) then + write (*, *) 'Lapack failed for forces' + stop + end if + forces_irred(:) = mB(1:n_irred, 1) + + ! Store the symmetrized forces + if (.not. allocated(vB)) allocate (vB(n)) + vB(:) = matmul(coeffM, forces_irred) + do i = 1, uc%na + v0 = vB((i - 1)*3 + 1:i*3) + v0 = matmul(uc%latticevectors, v0) + forces_sym(:, i) = lo_chop(v0, 1E-12_r8) + end do + end block fixforces + else + ! In case there are no internal degrees of freedom. + forces_sym = 0.0_r8 + end if +end subroutine + +!> expand a 3x3 symmetry operation to a 9x9 operation, to operate on flattened tensors +pure function expandoperation_pair(o) result(bigo) + !> the original 3x3 operation + real(r8), dimension(3, 3), intent(in) :: o + !> the resulting 9x9 operation + real(r8), dimension(9, 9) :: bigo + ! + integer :: i, j, ii, jj, k, l + ! + k = 0 + do i = 1, 3 + do ii = 1, 3 + k = k + 1 + l = 0 + do j = 1, 3 + do jj = 1, 3 + l = l + 1 + bigo(k, l) = o(ii, jj)*o(i, j) + end do + end do + end do + end do + bigo = lo_chop(bigo, 1E-11_r8) +end function + +end module diff --git a/src/libolle/Makefile b/src/libolle/Makefile index a1d75b1c..9110a283 100644 --- a/src/libolle/Makefile +++ b/src/libolle/Makefile @@ -563,7 +563,8 @@ $(OBJECT_PATH)type_forceconstant_thirdorder.o:\ $(OBJECT_PATH)gottochblandat.o\ $(OBJECT_PATH)type_crystalstructure.o\ $(OBJECT_PATH)type_qpointmesh.o\ - $(OBJECT_PATH)type_distancetable.o + $(OBJECT_PATH)type_distancetable.o\ + $(OBJECT_PATH)type_symmetryoperation.o $(FC) $(FFLAGS) -c type_forceconstant_thirdorder.f90 -o $@ $(OBJECT_PATH)type_forceconstant_fourthorder.o:\ $(OBJECT_PATH)konstanter.o\ @@ -639,6 +640,8 @@ $(OBJECT_PATH)type_embeddedatom.o:\ $(FC) $(FFLAGS) -c type_embeddedatom.f90 -o $@ $(OBJECT_PATH)type_phonon_dispersions.o:\ + type_phonon_dispersions.f90\ + type_phonon_dispersions_generation.f90\ $(OBJECT_PATH)konstanter.o\ $(OBJECT_PATH)gottochblandat.o\ $(OBJECT_PATH)mpi_wrappers.o\ diff --git a/src/libolle/ifc_solvers_lsq.f90 b/src/libolle/ifc_solvers_lsq.f90 index 962c0174..cf09be59 100644 --- a/src/libolle/ifc_solvers_lsq.f90 +++ b/src/libolle/ifc_solvers_lsq.f90 @@ -79,6 +79,54 @@ module subroutine lsq_solve(map, uc, ss, ih, tp, force_stable, usereference, max end if end block init + if (map%have_fc_singlet) then + first: block + type(lo_scaling_and_product) :: asc + real(r8), dimension(:, :), allocatable :: wA + real(r8), dimension(:), allocatable :: wB, wC + integer :: nrow + + ! Normal least squares for the full solution + call asc%generate(ih%cm1, ih%f1, nf, map%xuc%nx_fc_singlet, mw, noscale=.true.) + if (mw%r .eq. solrnk) then + call lo_linear_least_squares( & + asc%ATA, asc%ATB, map%xuc%x_fc_singlet, map%constraints%eq1, map%constraints%neq1, gramified=.true.) + end if + call mw%bcast(map%xuc%x_fc_singlet, solrnk, __FILE__, __LINE__) + call asc%destroy() + + ! Then the cross-checking guys + do div = 1, n_division + call tp%collect(map, 1, div, .true., nrow, mem, ih=ih, coeff=wA, values=wB) + call asc%generate(wA, wB, nrow, map%xuc%nx_fc_singlet, mw, noscale=.true.) + if (mw%r .eq. solrnk) then + call lo_linear_least_squares( & + asc%ATA, asc%ATB, ih%dx1(:, div), map%constraints%eq1, map%constraints%neq1, gramified=.true.) + end if + call asc%destroy() + call mem%deallocate(wA, persistent=.true., scalable=.true., file=__FILE__, line=__LINE__) + call mem%deallocate(wB, persistent=.true., scalable=.true., file=__FILE__, line=__LINE__) + end do + call mw%bcast(ih%dx1, solrnk, __FILE__, __LINE__) + + ! And subtract forces from the solution + if (tp%nt .gt. 0) then + call mem%allocate(wC, nf, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) + wC = 0.0_r8 + call lo_gemv(ih%cm1, map%xuc%x_fc_singlet, wC) + ih%f1 = wC + ih%f2 = ih%f2 - wC + ih%f3 = ih%f3 - wC + ih%f4 = ih%f4 - wC + call mem%deallocate(wC, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) + end if + if (verbosity .gt. 0) write (*, *) '... solved first order (', tochar(walltime() - timer), 's)' + end block first + else + ! No first order, set it to nothing. + ih%f1 = 0.0_r8 + end if + if (map%have_fc_pair) then second: block type(lo_scaling_and_product) :: asc @@ -131,7 +179,7 @@ module subroutine lsq_solve(map, uc, ss, ih, tp, force_stable, usereference, max call mem%allocate(wC, nf, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) wC = 0.0_r8 call lo_gemv(ih%cm2, map%xuc%x_fc_pair, wC) - ih%f1 = ih%f1 - wC + ! ih%f1 = ih%f1 - wC ih%f2 = wC ih%f3 = ih%f3 - wC ih%f4 = ih%f4 - wC @@ -148,54 +196,6 @@ module subroutine lsq_solve(map, uc, ss, ih, tp, force_stable, usereference, max ih%f2 = 0.0_r8 end if - if (map%have_fc_singlet) then - first: block - type(lo_scaling_and_product) :: asc - real(r8), dimension(:, :), allocatable :: wA - real(r8), dimension(:), allocatable :: wB, wC - integer :: nrow - - ! Normal least squares for the full solution - call asc%generate(ih%cm1, ih%f1, nf, map%xuc%nx_fc_singlet, mw, noscale=.true.) - if (mw%r .eq. solrnk) then - call lo_linear_least_squares( & - asc%ATA, asc%ATB, map%xuc%x_fc_singlet, map%constraints%eq1, map%constraints%neq1, gramified=.true.) - end if - call mw%bcast(map%xuc%x_fc_singlet, solrnk, __FILE__, __LINE__) - call asc%destroy() - - ! Then the cross-checking guys - do div = 1, n_division - call tp%collect(map, 1, div, .true., nrow, mem, ih=ih, coeff=wA, values=wB) - call asc%generate(wA, wB, nrow, map%xuc%nx_fc_singlet, mw, noscale=.true.) - if (mw%r .eq. solrnk) then - call lo_linear_least_squares( & - asc%ATA, asc%ATB, ih%dx1(:, div), map%constraints%eq1, map%constraints%neq1, gramified=.true.) - end if - call asc%destroy() - call mem%deallocate(wA, persistent=.true., scalable=.true., file=__FILE__, line=__LINE__) - call mem%deallocate(wB, persistent=.true., scalable=.true., file=__FILE__, line=__LINE__) - end do - call mw%bcast(ih%dx1, solrnk, __FILE__, __LINE__) - - ! And subtract forces from the solution - if (tp%nt .gt. 0) then - call mem%allocate(wC, nf, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) - wC = 0.0_r8 - call lo_gemv(ih%cm1, map%xuc%x_fc_singlet, wC) - ih%f1 = wC - !ih%f2=ih%f2-wC - ih%f3 = ih%f3 - wC - ih%f4 = ih%f4 - wC - call mem%deallocate(wC, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) - end if - if (verbosity .gt. 0) write (*, *) '... solved first order (', tochar(walltime() - timer), 's)' - end block first - else - ! No first order, set it to nothing. - ih%f1 = 0.0_r8 - end if - if (map%have_fc_triplet) then third: block type(lo_scaling_and_product) :: asc @@ -231,6 +231,7 @@ module subroutine lsq_solve(map, uc, ss, ih, tp, force_stable, usereference, max call mem%allocate(wC, nf, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) wC = 0.0_r8 call lo_gemv(ih%cm3, map%xuc%x_fc_triplet, wC) + ! ih%f1 = ih%f1 - wC ih%f3 = wC ih%f4 = ih%f4 - wC call mem%deallocate(wC, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) @@ -282,6 +283,7 @@ module subroutine lsq_solve(map, uc, ss, ih, tp, force_stable, usereference, max call mem%allocate(wC, nf, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) wC = 0.0_r8 call lo_gemv(ih%cm4, map%xuc%x_fc_quartet, wC) + ! ih%f1 = ih%f1 - wC ih%f4 = wC call mem%deallocate(wC, persistent=.false., scalable=.true., file=__FILE__, line=__LINE__) end if @@ -295,6 +297,7 @@ module subroutine lsq_solve(map, uc, ss, ih, tp, force_stable, usereference, max ! No fourth order, set it to nothing. ih%f4 = 0.0_r8 end if + end subroutine ! below are not exposed diff --git a/src/libolle/type_forceconstant_secondorder.f90 b/src/libolle/type_forceconstant_secondorder.f90 index 4ed7c70a..0e72ebde 100644 --- a/src/libolle/type_forceconstant_secondorder.f90 +++ b/src/libolle/type_forceconstant_secondorder.f90 @@ -6,7 +6,7 @@ module type_forceconstant_secondorder !! use konstanter, only: r8, i8, lo_iou, lo_tol, lo_sqtol, lo_pi, lo_twopi, lo_huge, lo_hugeint, lo_exitcode_param, lo_exitcode_symmetry, & lo_freqtol, lo_kb_hartree, lo_frequency_Hartree_to_THz -use gottochblandat, only: tochar, walltime, lo_clean_fractional_coordinates, lo_chop, lo_sqnorm +use gottochblandat, only: tochar, walltime, lo_clean_fractional_coordinates, lo_chop, lo_sqnorm, lo_invert_real_matrix use mpi_wrappers, only: lo_mpi_helper, lo_stop_gracefully use lo_memtracker, only: lo_mem_helper use lo_longrange_electrostatics, only: lo_ewald_parameters @@ -117,6 +117,7 @@ module type_forceconstant_secondorder !> elastic constants, could be fun or something real(r8), dimension(6, 6) :: elastic_constants_voigt = lo_huge real(r8), dimension(3, 3, 3, 3) :: elastic_constants_tensor = lo_huge + real(r8), dimension(6, 6) :: elastic_compliances_voigt = lo_huge ! Can be useful with the commensurate modes: @@ -373,6 +374,7 @@ subroutine remap(fc, uc, ss, fcss, ind) fcss%cutoff = -1 fcss%elastic_constants_voigt = -1 fcss%elastic_constants_tensor = -1 + fcss%elastic_compliances_voigt = -1 fcss%npairshells = -1 fcss%npairop = -1 fcss%nifc = -1 @@ -459,7 +461,7 @@ subroutine get_elastic_constants(fc, uc) integer :: i, j integer :: a1, pair integer :: al, be, gm, la - real(r8), dimension(6, 6) :: elc2 + real(r8), dimension(6, 6) :: elc2, s_ij real(r8), dimension(3, 3, 3, 3) :: elc, bracket real(r8), dimension(3, 3) :: m real(r8), dimension(3) :: r @@ -507,6 +509,12 @@ subroutine get_elastic_constants(fc, uc) fc%elastic_constants_voigt = elc2 fc%elastic_constants_tensor = elc + + ! add compliances (i.e. the inverse) + s_ij = 0.0_r8 + call lo_invert_real_matrix(elc2, s_ij) + fc%elastic_compliances_voigt = s_ij + contains !> Voigt-notation function contract_elastic_constant_indices(mu, nu) result(ind) diff --git a/src/libolle/type_forceconstant_thirdorder.f90 b/src/libolle/type_forceconstant_thirdorder.f90 index b7c8d70a..9057f6e7 100644 --- a/src/libolle/type_forceconstant_thirdorder.f90 +++ b/src/libolle/type_forceconstant_thirdorder.f90 @@ -6,6 +6,7 @@ module type_forceconstant_thirdorder use type_crystalstructure, only: lo_crystalstructure use type_distancetable, only: lo_distancetable use type_qpointmesh, only: lo_qpoint +use type_symmetryoperation, only: lo_operate_on_secondorder_tensor implicit none private @@ -60,6 +61,8 @@ module type_forceconstant_thirdorder procedure :: unallocate !> mode gruneisen parameter procedure :: mode_gruneisen_parameter + !> mode gruneisen tensor + procedure :: mode_gruneisen_tensor !> size in memory, in bytes procedure :: size_in_mem => fc3_size_in_mem end type @@ -500,14 +503,53 @@ subroutine mode_gruneisen_parameter(fct, p, qpoint, omega, eigenvectors, g) complex(r8), dimension(:, :), intent(in) :: eigenvectors !> the mode gruneisen parameters real(r8), dimension(:), intent(out) :: g + + integer ::nb, s + real(r8), dimension(:, :, :), allocatable :: g_tensor + + nb = fct%na*3 + allocate (g_tensor(nb, 3, 3)) + + ! get the tensors first + call mode_gruneisen_tensor(fct, p, qpoint, omega, eigenvectors, g_tensor) + + ! trace over tensor + do s = 1, nb + g(s) = (g_tensor(s, 1, 1) + g_tensor(s, 2, 2) + g_tensor(s, 3, 3))/3.0_r8 + ! g(s) = sum(g_tensor(s, :, :)) + end do + + deallocate (g_tensor) + +end subroutine + +!> Get the mode Gruneisen tensor from the third order force constants +subroutine mode_gruneisen_tensor(fct, p, qpoint, omega, eigenvectors, g_tensor) + !> the third order force constants + class(lo_forceconstant_thirdorder), intent(in) :: fct + !> The crystal structure + type(lo_crystalstructure), intent(in) :: p + !> the q-point + class(lo_qpoint), intent(in) :: qpoint + !> the frequencies + real(r8), dimension(:), intent(in) :: omega + !> the eigenvectors + complex(r8), dimension(:, :), intent(in) :: eigenvectors + !> the mode gruneisen tensors (mode,xyz,xyz) + real(r8), dimension(:, :, :), intent(out) :: g_tensor ! real(r8), dimension(3) :: qv, rv2, rv3 real(r8), dimension(3, 3, 3) :: m real(r8) :: effmass, iqr complex(r8), dimension(:, :, :), allocatable :: egv - complex(r8) :: c0, evprod, expiqr - integer :: i, j, k, s, t + complex(r8) :: evprod, expiqr + complex(r8), dimension(3, 3) :: c0_tensor + real(r8), dimension(3, 3) :: c0_tensor_real, g_ab_temp + integer :: i, j, k, s, t, l integer :: nb, atom1, atom2, atom3 + integer :: io, n_invariant_operation + real(r8), dimension(3) :: v0, v1, w0, w1 + complex(r8), dimension(3) :: cv0, cv1, cv2 nb = fct%na*3 qv = qpoint%r*lo_twopi @@ -525,11 +567,12 @@ subroutine mode_gruneisen_parameter(fct, p, qpoint, omega, eigenvectors, g) end do end do - g = 0.0_r8 + g_tensor = 0.0_r8 ! Get the actual gruneisen parameters do s = 1, nb ! A parameter for each band - c0 = 0.0_r8 + c0_tensor = 0.0_r8 + c0_tensor_real = 0.0_r8 do atom1 = 1, fct%na do t = 1, fct%atom(atom1)%n atom2 = fct%atom(atom1)%triplet(t)%i2 @@ -542,27 +585,56 @@ subroutine mode_gruneisen_parameter(fct, p, qpoint, omega, eigenvectors, g) expiqr = cmplx(cos(iqr), sin(iqr), r8) m = fct%atom(atom1)%triplet(t)%m do i = 1, 3 - do j = 1, 3 - evprod = conjg(egv(i, atom1, s))*egv(j, atom2, s) - do k = 1, 3 - c0 = c0 + m(i, j, k)*evprod*expiqr*rv3(k)*effmass + do j = 1, 3 + evprod = conjg(egv(i, atom1, s))*egv(j, atom2, s) + do k = 1, 3 + ! c0 = c0 + m(i, j, k)*evprod*expiqr*rv3(k)*effmass + do l = 1, 3 + c0_tensor(k, l) = c0_tensor(k, l) & + + 0.5_r8*evprod*expiqr*effmass & + *(m(i, j, k)*rv3(l) + m(i, j, l)*rv3(k)) + end do + end do end do end do - end do end do end do ! I don't want infinity at gamma if (omega(s) .lt. lo_freqtol) then - c0 = 0.0_r8 + c0_tensor = 0.0_r8 else - c0 = c0/(6.0_r8*(omega(s)**2)) + c0_tensor = c0_tensor/(2.0_r8*(omega(s)**2)) end if #ifdef AGRESSIVE_SANITY - if (abs(aimag(c0)) .gt. lo_sqtol) then + if (abs(aimag(sum(c0_tensor))) .gt. lo_tol) then + write (*, *) 'size: ', abs(aimag(sum(c0_tensor))) call lo_stop_gracefully(['Imaginary component of gruneisen parameters nonzero, not correct'], lo_exitcode_symmetry, __FILE__, __LINE__) end if #endif - g(s) = -real(c0) + + c0_tensor_real = -real(c0_tensor) + + ! FK: average over small group because we can + + g_ab_temp = 0.0_r8 + n_invariant_operation = qpoint%n_invariant_operation + if (n_invariant_operation >= 1) then + + do k = 1, qpoint%n_invariant_operation + io = qpoint%invariant_operation(k) + associate (op => p%sym%op(abs(io))) + g_ab_temp = g_ab_temp + lo_operate_on_secondorder_tensor(op, c0_tensor_real) + end associate + end do + g_ab_temp = g_ab_temp/real(n_invariant_operation, r8) + + else + ! nothing to do in this case + g_ab_temp = c0_tensor_real + + end if + + g_tensor(s, :, :) = g_ab_temp end do deallocate (egv) end subroutine diff --git a/src/libolle/type_mdsim.f90 b/src/libolle/type_mdsim.f90 index ecfb9662..98e6b9f8 100644 --- a/src/libolle/type_mdsim.f90 +++ b/src/libolle/type_mdsim.f90 @@ -7,7 +7,7 @@ module type_mdsim lo_time_fs_to_au, lo_time_au_to_fs use gottochblandat, only: open_file, lo_clean_fractional_coordinates, lo_mean, tochar, lo_trueNtimes, walltime, lo_sqnorm, & lo_progressbar_init, lo_progressbar, lo_find_rotation_that_makes_strain_diagonal, lo_trace, & - lo_mass_from_Z + lo_mass_from_Z, lo_stddev use mpi_wrappers, only: lo_mpi_helper, lo_stop_gracefully use hdf5_wrappers, only: lo_hdf5_helper, lo_h5_store_attribute, lo_h5_store_data, lo_h5_read_attribute, & lo_h5_read_data, lo_h5_does_dataset_exist @@ -1287,36 +1287,50 @@ subroutine read_from_file(sim, verbosity, stride, magnetic, dielectric, variable ! Calculate some things finalize: block + real(r8) :: n_samples real(r8), dimension(3, 3) :: m0 - real(r8) :: avg_pressure, avg_temperature + real(r8) :: pressure_avg, pressure_std, pressure_err, temperature_avg + real(r8), dimension(3, 3) :: stress_avg, stress_err integer :: i, j character(len=1000) :: opf - ! get averages of stat things - avg_pressure = lo_mean(sim%stat%pressure)*lo_pressure_HartreeBohr_to_GPa - avg_temperature = lo_mean(sim%stat%temperature) + + ! number of samples as float + n_samples = 1.0_r8*sim%nt + + ! get statistics including error + pressure_avg = lo_mean(-sim%stat%pressure)*lo_pressure_HartreeBohr_to_GPa + pressure_std = lo_stddev(-sim%stat%pressure)*lo_pressure_HartreeBohr_to_GPa + pressure_err = pressure_std/sqrt(n_samples) + temperature_avg = lo_mean(sim%stat%temperature) do i = 1, 3 do j = 1, 3 - m0(j, i) = lo_mean(sim%stat%stress(j, i, :)) + associate (stress => sim%stat%stress(j, i, :)) + stress_avg(j, i) = lo_mean(stress) + stress_err(j, i) = lo_stddev(stress)/sqrt((n_samples)) + end associate end do end do - m0 = m0*lo_pressure_HartreeBohr_to_GPa + stress_avg = stress_avg*lo_pressure_HartreeBohr_to_GPa + stress_err = stress_err*lo_pressure_HartreeBohr_to_GPa ! all done, perhaps dump some info if (verbosity .gt. 0) then - write (*, *) '... short summary of simulation:' + write (*, *) '... statistics for DFT pressure with 1.96x standard error (95% confidence interval):' write (*, *) ' number of atoms: ', tochar(sim%na) write (*, *) ' number of configurations: ', tochar(sim%nt) - write (*, *) ' average temperature (K): ', tochar(avg_temperature) + write (*, *) ' sqrt of number of configurations: ', tochar(sqrt(n_samples)) + write (*, *) ' average temperature (K): ', tochar(temperature_avg) write (*, *) ' thermostat set to (K): ', tochar(sim%temperature_thermostat) - write (*, *) ' average pressure (GPa): ', tochar(avg_pressure) - opf = "(1X,' average stresstensor xx xy xz (GPa): ',3(F12.4,' ') )" - write (*, opf) m0(:, 1) - opf = "(1X,' average stresstensor yx yy yz (GPa): ',3(F12.4,' ') )" - write (*, opf) m0(:, 2) - opf = "(1X,' average stresstensor zx zy zz (GPa): ',3(F12.4,' ') )" - write (*, opf) m0(:, 3) + write (*, *) ' average pressure (GPa): ', tochar(pressure_avg), ' +/- ', tochar(1.96*pressure_err) + opf = "(1X,' average stresstensor xx xy xz (GPa): ',3(F12.5,' '),' +/- ',3(F12.5,' ') )" + write (*, opf) stress_avg(:, 1), 1.96*stress_err(:, 1) + opf = "(1X,' average stresstensor yx yy yz (GPa): ',3(F12.5,' '),' +/- ',3(F12.5,' ') )" + write (*, opf) stress_avg(:, 2), 1.96*stress_err(:, 2) + opf = "(1X,' average stresstensor zx zy zz (GPa): ',3(F12.5,' '),' +/- ',3(F12.5,' ') )" + write (*, opf) stress_avg(:, 3), 1.96*stress_err(:, 3) write (*, *) '... finished reading simulation (', tochar(walltime() - timer), 's)' end if + end block finalize end subroutine diff --git a/src/libolle/type_phonon_dispersions.f90 b/src/libolle/type_phonon_dispersions.f90 index 5c759325..c8537ee6 100644 --- a/src/libolle/type_phonon_dispersions.f90 +++ b/src/libolle/type_phonon_dispersions.f90 @@ -6,7 +6,7 @@ module type_phonon_dispersions lo_sqtol, lo_temperaturetol, & lo_status, lo_hbar_hartree, lo_frequency_Hartree_to_Hz, lo_groupvel_HartreeBohr_to_ms, & lo_exitcode_param, lo_hartree_to_eV, lo_bohr_to_m, lo_kappa_au_to_SI, lo_kb_hartree, & - lo_Hartree_to_Joule, lo_exitcode_physical, lo_imag, lo_time_s_to_au + lo_Hartree_to_Joule, lo_exitcode_physical, lo_imag, lo_time_s_to_au, lo_volume_bohr_to_A use gottochblandat, only: open_file, walltime, tochar, lo_trueNtimes, lo_classical_harmonic_oscillator_free_energy, & lo_harmonic_oscillator_cv, lo_harmonic_oscillator_entropy, lo_harmonic_oscillator_free_energy, & lo_planck, lo_sqnorm, lo_progressbar_init, lo_progressbar, qsort, lo_outerproduct, lo_chop, & @@ -21,7 +21,8 @@ module type_phonon_dispersions use type_forceconstant_thirdorder, only: lo_forceconstant_thirdorder use type_qpointmesh, only: lo_qpoint, lo_qpoint_mesh, lo_monkhorst_pack_mesh, lo_wedge_mesh, lo_fft_mesh, & lo_get_small_group_of_qpoint -use type_symmetryoperation, only: lo_eigenvector_transformation_matrix, lo_operate_on_vector, lo_expandoperation_pair +use type_symmetryoperation, only: lo_eigenvector_transformation_matrix, lo_operate_on_vector, lo_expandoperation_pair, & + lo_operate_on_secondorder_tensor use hdf5_wrappers, only: lo_hdf5_helper implicit none @@ -37,8 +38,9 @@ module type_phonon_dispersions real(r8), dimension(:, :), allocatable :: vel !> mode eigenvectors complex(r8), dimension(:, :), allocatable :: egv - !> mode gruneisen parameter + !> mode gruneisen parameter and tensor real(r8), dimension(:), allocatable :: gruneisen + real(r8), dimension(:, :, :), allocatable :: gruneisen_tensor !> linewidth real(r8), dimension(:), allocatable :: linewidth !> anharmonic frequency shift @@ -125,7 +127,8 @@ module type_phonon_dispersions !> create the full dispersions procedure :: generate !> calculate gruneisen parameters - !procedure :: gruneisen + procedure :: gruneisen + procedure :: gruneisen_tensor !> phonon free energy procedure :: phonon_free_energy !> classical phonon free energy @@ -163,14 +166,22 @@ module subroutine generate(dr, qp, fc, p, mw, mem, verbosity) type(lo_mem_helper), intent(inout) :: mem integer, intent(in) :: verbosity end subroutine - ! module subroutine gruneisen(dr,qp,fct,p,verbosity,mpi_communicator) - ! class(lo_phonon_dispersions), intent(inout) :: dr - ! class(lo_qpoint_mesh), intent(inout) :: qp - ! type(lo_forceconstant_thirdorder), intent(in) :: fct - ! type(lo_crystalstructure), intent(in) :: p - ! integer, intent(in), optional :: verbosity - ! integer, intent(in), optional :: mpi_communicator - ! end subroutine + module subroutine gruneisen(dr, qp, fct, p, verbosity, mpi_communicator) + class(lo_phonon_dispersions), intent(inout) :: dr + class(lo_qpoint_mesh), intent(inout) :: qp + type(lo_forceconstant_thirdorder), intent(in) :: fct + type(lo_crystalstructure), intent(in) :: p + integer, intent(in), optional :: verbosity + integer, intent(in), optional :: mpi_communicator + end subroutine + module subroutine gruneisen_tensor(dr, qp, fct, p, verbosity, mpi_communicator) + class(lo_phonon_dispersions), intent(inout) :: dr + class(lo_qpoint_mesh), intent(inout) :: qp + type(lo_forceconstant_thirdorder), intent(in) :: fct + type(lo_crystalstructure), intent(in) :: p + integer, intent(in), optional :: verbosity + integer, intent(in), optional :: mpi_communicator + end subroutine module subroutine allgather_irreducible(dr, mw, mem) class(lo_phonon_dispersions), intent(inout) :: dr type(lo_mpi_helper), intent(inout) :: mw @@ -331,12 +342,13 @@ subroutine write_to_hdf5(dr, qp, uc, filename, mem, temperature) real(r8), dimension(:, :), allocatable :: dd real(r8), dimension(3) :: v0, v1 real(r8) :: n1, n, f0, omega, cv, tau - integer :: i, j, k, l + integer :: i, j, k, l, s character(len=1000) :: dname call h5%init(__FILE__, __LINE__) call h5%open_file('write', trim(filename)) ! Some general attributes first: + call h5%store_attribute(uc%volume*lo_volume_bohr_to_A, h5%file_id, 'volume', lo_status) call h5%store_attribute(dr%n_mode/3, h5%file_id, 'number_of_atoms', lo_status) call h5%store_attribute(dr%n_mode, h5%file_id, 'number_of_bands', lo_status) call h5%store_attribute(dr%n_full_qpoint, h5%file_id, 'number_of_qpoints', lo_status) @@ -408,6 +420,19 @@ subroutine write_to_hdf5(dr, qp, uc, filename, mem, temperature) call h5%store_data(dd, h5%file_id, trim(dname), enhet='dimensionless', dimensions='q-vector,mode') call mem%deallocate(dd, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) end if + if (allocated(dr%aq(1)%gruneisen_tensor)) then + dname = 'gruneisen_tensors' + call mem%allocate(dddd, [3, 3, dr%n_mode, qp%n_full_point], & + persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + dddd = 0.0_r8 + do i = 1, qp%n_full_point + do s = 1, dr%n_mode + dddd(:, :, s, i) = dr%aq(i)%gruneisen_tensor(s, :, :) + end do + end do + call h5%store_data(dddd, h5%file_id, trim(dname), enhet='dimensionless', dimensions='q-vector,mode,xyz,xyz') + call mem%deallocate(dddd, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) + end if ! Lifetimes if (allocated(dr%iq(1)%linewidth)) then @@ -1192,6 +1217,7 @@ function phonon_dispersions_qpoint_size_in_mem(p) result(mem) if (allocated(p%vel)) mem = mem + storage_size(p%vel)*size(p%vel) if (allocated(p%egv)) mem = mem + storage_size(p%egv)*size(p%egv) if (allocated(p%gruneisen)) mem = mem + storage_size(p%gruneisen)*size(p%gruneisen) + if (allocated(p%gruneisen_tensor)) mem = mem + storage_size(p%gruneisen_tensor)*size(p%gruneisen_tensor) if (allocated(p%linewidth)) mem = mem + storage_size(p%linewidth)*size(p%linewidth) if (allocated(p%shift3)) mem = mem + storage_size(p%shift3)*size(p%shift3) if (allocated(p%shift4)) mem = mem + storage_size(p%shift4)*size(p%shift4) @@ -1230,6 +1256,7 @@ pure function phonon_dispersions_qpoint_size_packed(p) result(mem) if (allocated(p%vel)) mem = mem + storage_size(p%vel)*size(p%vel)/8 if (allocated(p%egv)) mem = mem + storage_size(p%egv)*size(p%egv)/8 if (allocated(p%gruneisen)) mem = mem + storage_size(p%gruneisen)*size(p%gruneisen)/8 + if (allocated(p%gruneisen_tensor)) mem = mem + storage_size(p%gruneisen_tensor)*size(p%gruneisen_tensor)/8 if (allocated(p%linewidth)) mem = mem + storage_size(p%linewidth)*size(p%linewidth)/8 if (allocated(p%shift3)) mem = mem + storage_size(p%shift3)*size(p%shift3)/8 if (allocated(p%shift4)) mem = mem + storage_size(p%shift4)*size(p%shift4)/8 diff --git a/src/libolle/type_phonon_dispersions_generation.f90 b/src/libolle/type_phonon_dispersions_generation.f90 index 556e9025..2cde45ee 100644 --- a/src/libolle/type_phonon_dispersions_generation.f90 +++ b/src/libolle/type_phonon_dispersions_generation.f90 @@ -358,101 +358,198 @@ module subroutine allgather_fullmesh(dr, mw, mem) call mem%deallocate(rbuf, persistent=.false., scalable=.false., file=__FILE__, line=__LINE__) end subroutine -! !> Calculate mode gruneisen parameters on a q-point mesh -! module subroutine gruneisen(dr,qp,fct,p,verbosity,mpi_communicator) -! !> the dispersion relations -! class(lo_phonon_dispersions), intent(inout) :: dr -! !> the qpoint mesh -! class(lo_qpoint_mesh), intent(inout) :: qp -! !> the third order force constant -! type(lo_forceconstant_thirdorder), intent(in) :: fct -! !> the crystal structure -! type(lo_crystalstructure), intent(in) :: p -! !> how much to talk -! integer, intent(in), optional :: verbosity -! !> do the calculation distributed -! integer, intent(in), optional :: mpi_communicator -! ! -! write(*,*) 'FIXME GRUNEISEN MPI' -! stop -! -! integer :: q,rank,nrank,lqp -! real(r8) :: t0 -! logical :: usempi -! -! t0=walltime() -! if ( present(verbosity) ) then -! verbosity=verbosity -! else -! verbosity=0 -! endif -! -! ! If we have an MPI communicator attached, check that the q-points are distributed accordingly. -! if ( present(mpi_communicator) ) then -! ! Are the q-points distributed across a communicator? -! if ( qp%mpi%initialized ) then -! ! Yep, distributed. The same communicator that is sent here? -! if ( qp%mpi%communicator .ne. mpi_communicator ) then -! ! This is weird. Not sure what to do. Probably redistribute? -! call qp%divide_across_mpi(mpi_communicator) -! endif -! else -! ! Nope, not initialized. Do that! -! call qp%divide_across_mpi(mpi_communicator) -! endif -! usempi=.true. -! call mpi_comm_rank(mpi_communicator,rank,lo_status) -! call mpi_comm_size(mpi_communicator,nrank,lo_status) -! ! remove verbosity except on head rank -! if ( rank .ne. 0 ) verbosity=0 -! ! if only one rank, don't use MPI at all. -! if ( nrank .eq. 1 ) usempi=.false. -! else -! ! Don't use MPI at all. -! usempi=.false. -! rank=0 -! nrank=1 -! endif -! -! if ( verbosity .gt. 0 ) then -! if ( usempi ) then -! write(*,*) '... calculating gruneisen parameters over ',tochar(nrank),' MPI ranks' -! else -! write(*,*) '... calculating gruneisen parameters' -! endif -! endif -! -! ! Make some space -! do q=1,qp%n_irr_point -! lo_allocate(dr%iq(q)%gruneisen(dr%n_mode)) -! dr%iq(q)%gruneisen=0.0_r8 -! enddo -! do q=1,qp%n_full_point -! lo_allocate(dr%aq(q)%gruneisen(dr%n_mode)) -! dr%aq(q)%gruneisen=0.0_r8 -! enddo -! -! ! And the actual calculation -! if ( usempi ) then -! do lqp=1,qp%mpi%nqirr -! q=qp%mpi%qirr(lqp) -! call fct%mode_gruneisen_parameter(p,qp%ip(q),dr%iq(q)%omega,dr%iq(q)%egv,dr%iq(q)%gruneisen) -! enddo -! call dr%allgather_irreducible(qp,mpi_communicator,fieldlist=[4]) -! else -! do q=1,qp%n_irr_point -! call fct%mode_gruneisen_parameter(p,qp%ip(q),dr%iq(q)%omega,dr%iq(q)%egv,dr%iq(q)%gruneisen) -! enddo -! endif -! -! ! and spread them to the full thing -! do q=1,qp%n_full_point -! dr%aq(q)%gruneisen=dr%iq( qp%ap(q)%irrind )%gruneisen -! enddo -! if ( verbosity .gt. 0 ) then -! write(*,*) '... got gruneisen parameters in irreducible wedge in '//tochar(walltime()-t0,6)//' s' -! t0=walltime() -! endif -! end subroutine +!> Calculate mode gruneisen parameters on a q-point mesh +module subroutine gruneisen(dr, qp, fct, p, verbosity, mpi_communicator) + !> the dispersion relations + class(lo_phonon_dispersions), intent(inout) :: dr + !> the qpoint mesh + class(lo_qpoint_mesh), intent(inout) :: qp + !> the third order force constant + type(lo_forceconstant_thirdorder), intent(in) :: fct + !> the crystal structure + type(lo_crystalstructure), intent(in) :: p + !> how much to talk + integer, intent(in), optional :: verbosity + integer :: talk + !> do the calculation distributed + integer, intent(in), optional :: mpi_communicator + !! fkdev: let's go without MPI + ! + ! write (*, *) 'FIXME GRUNEISEN MPI' + ! stop + + integer :: q, rank, nrank, lqp + real(r8) :: t0 + logical :: usempi + + ! fkev: MPI currently not supported + if (present(mpi_communicator)) then + write (*, *) 'FIXME GRUNEISEN MPI' + stop + end if + + t0 = walltime() + if (present(verbosity)) then + talk = verbosity + else + talk = 0 + end if + + !! fkdev: comment out everything MPI related + ! ! If we have an MPI communicator attached, check that the q-points are distributed accordingly. + ! if (present(mpi_communicator)) then + ! ! Are the q-points distributed across a communicator? + ! if (qp%mpi%initialized) then + ! ! Yep, distributed. The same communicator that is sent here? + ! if (qp%mpi%communicator .ne. mpi_communicator) then + ! ! This is weird. Not sure what to do. Probably redistribute? + ! call qp%divide_across_mpi(mpi_communicator) + ! end if + ! else + ! ! Nope, not initialized. Do that! + ! call qp%divide_across_mpi(mpi_communicator) + ! end if + ! usempi = .true. + ! call mpi_comm_rank(mpi_communicator, rank, lo_status) + ! call mpi_comm_size(mpi_communicator, nrank, lo_status) + ! ! remove verbosity except on head rank + ! if (rank .ne. 0) _verbosity = 0 + ! ! if only one rank, don't use MPI at all. + ! if (nrank .eq. 1) usempi = .false. + ! else + ! ! Don't use MPI at all. + ! usempi = .false. + ! rank = 0 + ! nrank = 1 + ! end if + + !! fkdev: no mpi + usempi = .false. + rank = 0 + nrank = 1 + + if (talk .gt. 0) then + if (usempi) then + write (*, *) '... calculating gruneisen parameters over ', tochar(nrank), ' MPI ranks' + else + write (*, *) '... calculating gruneisen parameters' + end if + end if + + ! Make some space + do q = 1, qp%n_irr_point + allocate (dr%iq(q)%gruneisen(dr%n_mode)) + dr%iq(q)%gruneisen = 0.0_r8 + end do + do q = 1, qp%n_full_point + allocate (dr%aq(q)%gruneisen(dr%n_mode)) + dr%aq(q)%gruneisen = 0.0_r8 + end do + + ! And the actual calculation + if (usempi) then + ! do lqp = 1, qp%mpi%nqirr + ! q = qp%mpi%qirr(lqp) + ! call fct%mode_gruneisen_parameter(p, qp%ip(q), dr%iq(q)%omega, dr%iq(q)%egv, dr%iq(q)%gruneisen) + ! end do + ! call dr%allgather_irreducible(qp, mpi_communicator, fieldlist=[4]) + else + do q = 1, qp%n_irr_point + call fct%mode_gruneisen_parameter(p, qp%ip(q), dr%iq(q)%omega, dr%iq(q)%egv, dr%iq(q)%gruneisen) + end do + end if + + ! and spread them to the full thing + do q = 1, qp%n_full_point + dr%aq(q)%gruneisen = dr%iq(qp%ap(q)%irreducible_index)%gruneisen + end do + if (talk .gt. 0) then + write (*, *) '... got gruneisen parameters in irreducible wedge in '//tochar(walltime() - t0, 6)//' s' + t0 = walltime() + end if +end subroutine + +!> Calculate mode gruneisen tensors on a q-point mesh +module subroutine gruneisen_tensor(dr, qp, fct, p, verbosity, mpi_communicator) + !> the dispersion relations + class(lo_phonon_dispersions), intent(inout) :: dr + !> the qpoint mesh + class(lo_qpoint_mesh), intent(inout) :: qp + !> the third order force constant + type(lo_forceconstant_thirdorder), intent(in) :: fct + !> the crystal structure + type(lo_crystalstructure), intent(in) :: p + !> how much to talk + integer, intent(in), optional :: verbosity + integer :: talk + !> do the calculation distributed + integer, intent(in), optional :: mpi_communicator + + integer :: q, rank, nrank, lqp, iqp, iop, s + real(r8) :: t0 + logical :: usempi + + ! fkev: MPI currently not supported + if (present(mpi_communicator)) then + write (*, *) 'FIXME GRUNEISEN MPI' + stop + end if + + t0 = walltime() + if (present(verbosity)) then + talk = verbosity + else + talk = 0 + end if + + !! fkdev: no mpi + usempi = .false. + rank = 0 + nrank = 1 + + if (talk .gt. 0) then + write (*, *) '... calculating gruneisen parameters' + end if + + ! Make some space + do q = 1, qp%n_irr_point + allocate (dr%iq(q)%gruneisen_tensor(dr%n_mode, 3, 3)) + dr%iq(q)%gruneisen_tensor = 0.0_r8 + end do + do q = 1, qp%n_full_point + allocate (dr%aq(q)%gruneisen_tensor(dr%n_mode, 3, 3)) + dr%aq(q)%gruneisen_tensor = 0.0_r8 + ! call fct%mode_gruneisen_tensor(p, qp%ap(q), dr%aq(q)%omega, dr%aq(q)%egv, dr%aq(q)%gruneisen_tensor) + end do + + ! FK: we want to use symmetry like a human + ! compute in the irreducible wedge + do q = 1, qp%n_irr_point + call fct%mode_gruneisen_tensor(p, qp%ip(q), dr%iq(q)%omega, dr%iq(q)%egv, dr%iq(q)%gruneisen_tensor) + end do + + ! and spread them to the full thing + do q = 1, qp%n_full_point + iqp = qp%ap(q)%irreducible_index + iop = qp%ap(q)%operation_from_irreducible + associate ( & + ig_ab => dr%iq(iqp)%gruneisen_tensor, & + op => p%sym%op(abs(iop)), & + nb => fct%na*3 & + ) + ! we have a 3x3 matrix for each mode s that we need to rotate + do s = 1, nb + dr%aq(q)%gruneisen_tensor(s, :, :) = & + lo_operate_on_secondorder_tensor(op, ig_ab(s, :, :)) + end do + end associate + + end do + + if (talk .gt. 0) then + write (*, *) '... got gruneisen tensors in irreducible wedge in '//tochar(walltime() - t0, 6)//' s' + t0 = walltime() + end if +end subroutine end submodule diff --git a/src/phonon_dispersion_relations/main.f90 b/src/phonon_dispersion_relations/main.f90 index 314860e6..3987141c 100644 --- a/src/phonon_dispersion_relations/main.f90 +++ b/src/phonon_dispersion_relations/main.f90 @@ -318,7 +318,9 @@ program phonon_dispersion_relations stop end if if (opts%gruneisen) then - !call dr%gruneisen(qp,fct,uc,opts%verbosity) + call dr%gruneisen(qp, fct, uc, opts%verbosity) + ! fkdev: this means we do the work twice, could be nicer + call dr%gruneisen_tensor(qp, fct, uc, opts%verbosity) end if ! actually write it if (mw%talk) then From 15b110e430f5ef97bf42c342a30e1314a6781144 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Mon, 2 Oct 2023 10:47:17 +0200 Subject: [PATCH 2/3] mdsim | compute pressure from stress, just to be sure --- src/libolle/type_mdsim.f90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/libolle/type_mdsim.f90 b/src/libolle/type_mdsim.f90 index 98e6b9f8..9959dc84 100644 --- a/src/libolle/type_mdsim.f90 +++ b/src/libolle/type_mdsim.f90 @@ -1128,6 +1128,9 @@ subroutine read_from_file(sim, verbosity, stride, magnetic, dielectric, variable sim%stat%stress(3, 2, i) = sigma(5) sim%stat%stress(1, 3, i) = sigma(6) sim%stat%stress(3, 1, i) = sigma(6) + ! we do not trust the input file and make sure that + ! p = - 1/3 Tr(stress) + sim%stat%pressure(i) = -lo_trace(sim%stat%stress(:, :, i))/3.0_r8 else read (u, *) end if @@ -1299,8 +1302,8 @@ subroutine read_from_file(sim, verbosity, stride, magnetic, dielectric, variable n_samples = 1.0_r8*sim%nt ! get statistics including error - pressure_avg = lo_mean(-sim%stat%pressure)*lo_pressure_HartreeBohr_to_GPa - pressure_std = lo_stddev(-sim%stat%pressure)*lo_pressure_HartreeBohr_to_GPa + pressure_avg = lo_mean(sim%stat%pressure)*lo_pressure_HartreeBohr_to_GPa + pressure_std = lo_stddev(sim%stat%pressure)*lo_pressure_HartreeBohr_to_GPa pressure_err = pressure_std/sqrt(n_samples) temperature_avg = lo_mean(sim%stat%temperature) do i = 1, 3 From bb6e8f994a87ecc870444d0435305aa36269e7d0 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Mon, 2 Oct 2023 10:47:53 +0200 Subject: [PATCH 3/3] pressure | make sure not to forget about polar contributions --- src/extract_forceconstants/main.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/extract_forceconstants/main.f90 b/src/extract_forceconstants/main.f90 index 497da97f..ea4013b4 100644 --- a/src/extract_forceconstants/main.f90 +++ b/src/extract_forceconstants/main.f90 @@ -1102,8 +1102,8 @@ program extract_forceconstants ! do everything else in GPa for digestive purposes s0t = s0t*lo_pressure_HartreeBohr_to_GPa - s2t = s2t*lo_pressure_HartreeBohr_to_GPa spt = spt*lo_pressure_HartreeBohr_to_GPa + s2t = s2t*lo_pressure_HartreeBohr_to_GPa s3t = s3t*lo_pressure_HartreeBohr_to_GPa s4t = s4t*lo_pressure_HartreeBohr_to_GPa srt = srt*lo_pressure_HartreeBohr_to_GPa @@ -1123,9 +1123,9 @@ program extract_forceconstants end do p0t = p0t - s0t(i, i, :)/3.0_r8 p2t = p2t - s2t(i, i, :)/3.0_r8 - ppt = ppt - spt(i, i, :)/3.0_r8 p3t = p3t - s3t(i, i, :)/3.0_r8 p4t = p4t - s4t(i, i, :)/3.0_r8 + ppt = ppt - spt(i, i, :)/3.0_r8 prt = prt - srt(i, i, :)/3.0_r8 end do @@ -1139,9 +1139,9 @@ program extract_forceconstants ! pressure p0 = lo_mean(p0t) p2 = lo_mean(p2t) - pp = lo_mean(ppt) p3 = lo_mean(p3t) p4 = lo_mean(p4t) + pp = lo_mean(ppt) pr = lo_mean(prt) p0_std = lo_stddev(p0t) p2_std = lo_stddev(p2t) @@ -1175,15 +1175,15 @@ program extract_forceconstants write (*, *) '' write (*, "(1X,A,11X,A,10X,A,7X,A)") 'ENERGIES:', 'rms', 'stddev', 'stddev(residual) (meV/atom)' write (*, '(1X,A15,2(1X,F12.6),5X,A)') ' input:', sqrt(lo_mean(e0t**2))*tomev, lo_stddev(e0t)*tomev, '-' - write (*, '(1X,A15,3(1X,F12.6))') 'second order:', sqrt(lo_mean(e2t**2))*tomev, lo_stddev(e2t)*tomev, lo_stddev(e0t - e2t)*tomev write (*, '(1X,A15,3(1X,F12.6))') ' polar:', sqrt(lo_mean(ept**2))*tomev, lo_stddev(ept)*tomev, lo_stddev(e0t - ept - e2t)*tomev + write (*, '(1X,A15,3(1X,F12.6))') 'second order:', sqrt(lo_mean(e2t**2))*tomev, lo_stddev(e2t)*tomev, lo_stddev(e0t - e2t)*tomev write (*, '(1X,A15,3(1X,F12.6))') 'third order:', sqrt(lo_mean(e3t**2))*tomev, lo_stddev(e3t)*tomev, lo_stddev(e0t - ept - e2t - e3t)*tomev write (*, '(1X,A15,3(1X,F12.6))') 'fourth order:', sqrt(lo_mean(e4t**2))*tomev, lo_stddev(e4t)*tomev, lo_stddev(e0t - ept - e2t - e3t - e4t)*tomev write (*, "(1X,A,11X,A,10X,A,5X,A)") 'PRESSURE: ', 'mean', 'stddev', 'stdev (residual) (GPa)' write (*, '(1X,A15,2(1X,F12.6),5X,A)') ' input:', p0, p0_std, '-' - write (*, '(1X,A15,3(1X,F12.6),5X,A)') '2nd order:', p2, p2_std, lo_stddev(p0t - p2t) write (*, '(1X,A15,3(1X,F12.6),5X,A)') ' polar:', pp, pp_std, lo_stddev(p0t - p2t - ppt) + write (*, '(1X,A15,3(1X,F12.6),5X,A)') '2nd order:', p2, p2_std, lo_stddev(p0t - p2t) write (*, '(1X,A15,3(1X,F12.6),5X,A)') '3rd order:', p3, p3_std, lo_stddev(p0t - p2t - ppt - p3t) write (*, '(1X,A15,3(1X,F12.6),5X,A)') '4th order:', p4, p4_std, lo_stddev(p0t - p2t - ppt - p3t - p4t) write (*, '(1X,A15,2(1X,F12.6),5X,A)') ' residual:', pr, pr_std, '-' @@ -1191,14 +1191,14 @@ program extract_forceconstants write (*, *) associate (prefactor => 2.0_r8/sim%na/lo_kb_Hartree/3.0_r8) write (*, "(1X,A,F12.6,A,F12.6,A)") 'Harmonic temperature: ', & - prefactor*lo_mean(e2t), ' +/- ', prefactor*lo_stddev(e2t)/sqrt(n_samples), ' K' + prefactor*lo_mean(e2t + ept), ' +/- ', prefactor*lo_stddev(e2t + ept)/sqrt(n_samples), ' K' end associate write (*, *) associate (prefactor => 2.0_r8/3.0_r8/volume*lo_pressure_HartreeBohr_to_GPa) write (*, "(1X,A)") 'Harmonic pressure from energy p = 2/3 E_pot / V: ' write (*, "(1X,F12.6,A,F12.6,A)") & - prefactor*lo_mean(e2t), ' +/- ', prefactor*lo_stddev(e2t)/sqrt(n_samples), ' GPa' + prefactor*lo_mean(e2t + ept), ' +/- ', prefactor*lo_stddev(e2t + ept)/sqrt(n_samples), ' GPa' end associate write (*, *)