From 1efe2a31cd719047180064b41cb52a55e29cda44 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Thu, 12 Dec 2024 13:41:43 +0100 Subject: [PATCH 01/11] average_structure | add from tdep-devel --- build_things.sh | 1 + src/average_structure/Makefile | 33 + src/average_structure/main.f90 | 245 +++++ src/average_structure/options.f90 | 69 ++ src/average_structure/type_irredcell.f90 | 1031 ++++++++++++++++++ src/average_structure/type_irredcell_aux.f90 | 430 ++++++++ 6 files changed, 1809 insertions(+) create mode 100644 src/average_structure/Makefile create mode 100644 src/average_structure/main.f90 create mode 100644 src/average_structure/options.f90 create mode 100644 src/average_structure/type_irredcell.f90 create mode 100644 src/average_structure/type_irredcell_aux.f90 diff --git a/build_things.sh b/build_things.sh index 0b6d4e01..19c71ebe 100755 --- a/build_things.sh +++ b/build_things.sh @@ -146,6 +146,7 @@ refine_structure thermal_conductivity thermal_conductivity_2023 anharmonic_free_energy +average_structure " # only when we have cgal diff --git a/src/average_structure/Makefile b/src/average_structure/Makefile new file mode 100644 index 00000000..df6149ee --- /dev/null +++ b/src/average_structure/Makefile @@ -0,0 +1,33 @@ +include Makefile.inc +CODE = average_structure +PROG = ../../build/$(CODE)/$(CODE) +OBJECT_PATH=../../build/$(CODE)/ + +OBJS = $(OBJECT_PATH)main.o $(OBJECT_PATH)options.o $(OBJECT_PATH)type_irredcell.o + +LPATH = -L../../lib $(blaslapackLPATH) $(incLPATHhdf) $(incLPATHmpi) +IPATH = -I../../inc/libolle -I../../inc/libflap $(blaslapackIPATH) $(incIPATHhdf) $(incIPATHmpi) +LIBS = -lolle -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) + +OPT = -O0 -fcheck=all -fbacktrace -fno-realloc-lhs -finit-real=nan -finit-derived -fmax-errors=10 +F90 = $(FC) $(LPATH) $(IPATH) $(MODULE_FLAG) $(OBJECT_PATH) +F90FLAGS = $(OPT) $(MODS) $(LIBS) + +all: $(PROG) + +$(PROG): $(OBJS) + $(F90) $(LDFLAGS) -o $@ $(OBJS) $(LIBS) + +clean: + rm -f $(PROG) $(OBJS) *.mod + +.f90.o: + $(F90) $(F90FLAGS) -c $< $(LIBS) + +$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o $(OBJECT_PATH)type_irredcell.o + $(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@ +$(OBJECT_PATH)type_irredcell.o: + $(F90) $(F90FLAGS) -c type_irredcell.f90 $(LIBS) -o $@ +$(OBJECT_PATH)options.o: + $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ + diff --git a/src/average_structure/main.f90 b/src/average_structure/main.f90 new file mode 100644 index 00000000..99d3acf4 --- /dev/null +++ b/src/average_structure/main.f90 @@ -0,0 +1,245 @@ +#include "precompilerdefinitions" +program average_structure +!!{!holding/average_structure/manual.md!} +use konstanter !, only: flyt,lo_status,lo_volume_A_to_Bohr,lo_volume_bohr_to_A,lo_eV_to_Hartree,lo_Hartree_to_eV,lo_pressure_HartreeBohr_to_GPa +use gottochblandat +use options, only: lo_opts +use type_crystalstructure, only: lo_crystalstructure +use type_irredcell !, only: lo_irredcell +use type_mdsim, only: lo_mdsim +use type_symmetryoperation, only: lo_symset, lo_operate_on_vector +implicit none +! for normal phonon things +type(lo_opts) :: opts +type(lo_crystalstructure) :: uc, ss, ur, nuc, nss +type(lo_irredcell) :: ic +type(lo_mdsim) :: sim +real(flyt), dimension(:), allocatable :: volume, energy, eta, pressure, bulkmodulus +real(flyt), dimension(:), allocatable :: energy_model +integer :: u, i, ne +character(len=1000) :: opf + +init: block + integer :: u + integer :: i, j, k + real(flyt), dimension(6) :: sigma + character(len=1000) :: opf + real(flyt), dimension(3, 3) :: m0 + ! Get the command line arguments + call opts%parse() + + ! Get the unitcell that determines the symmetry + call uc%readfromfile('infile.ucposcar') + ! Get the unitcell that defines reference positions + call ss%readfromfile('infile.ssposcar') + call ss%classify('supercell', uc) + + ! Reference positions to determine displacement from + call ur%readfromfile('infile.refposcar') + + ! Read the MD simulation + call sim%read_from_file(verbosity=2, stride=opts%stride) +! lo_deallocate(sim%lattice) + ! Remove center of mass drift +! call sim%remove_force_and_center_of_mass_drift() + + call ic%generate(uc, verbosity=2, reference_positions=ur%r) + call avg_structure_nvt(uc, ss, ic, sim, nuc, nss, ur) + +!stop +! call sim%remove_force_and_center_of_mass_drift() +! ! Remove arbitrary rotations +! nuc=uc +! nss=ss +! do i=1,1 !20 +! m0=nuc%latticevectors +! ! Get the symmetric representation +! if ( i .eq. 1 ) then +! call ic%generate(nuc,verbosity=2) +! else +! call ic%generate(nuc,verbosity=0) +! endif +! ! +!! call remove_rotations(nuc,nss,ic,sim) +! ! Get the averaged structure +! uc=nuc +! ss=nss +! call avg_structure(uc,ss,ic,sim,nuc,nss) +! ! +! write(*,*) i,sum(abs(m0-nuc%latticevectors)) +! if ( sum(abs(m0-nuc%latticevectors)) .lt. 1E-10_flyt ) exit +! enddo + + ! Dump the new, clean averaged structure + call nuc%writetofile('outfile.ucposcar', 1) + call nss%writetofile('outfile.ssposcar', 1) + + write (*, *) '... wrote averaged structure' + + stop + +! u=open_file('out','outfile.meta') +! write(u,*) nss%na +! write(u,*) sim%nt +! write(u,*) sim%ts +! write(u,*) sim%t_thermostat +! close(u) +! write(*,*) '... wrote new stat' +! +! u=open_file('out','outfile.positions') +! do i=1,sim%nt +! do j=1,sim%na +! write(u,*) sim%r(:,j,i) +! enddo +! enddo +! close(u) +! write(*,*) '... wrote new positions' +! +! u=open_file('out','outfile.forces') +! do i=1,sim%nt +! do j=1,sim%na +! write(u,*) sim%f(:,j,i)*lo_force_HartreeBohr_to_eVA +! enddo +! enddo +! close(u) +! write(*,*) '... wrote new forces' +! +! u=open_file('out','outfile.stat') +! do i=1,sim%nt +! sigma(1)=sim%stat%stress(1,1,i) +! sigma(2)=sim%stat%stress(2,2,i) +! sigma(3)=sim%stat%stress(3,3,i) +! sigma(4)=sim%stat%stress(1,2,i) +! sigma(5)=sim%stat%stress(2,3,i) +! sigma(6)=sim%stat%stress(1,3,i) +! sigma=sigma*lo_pressure_HartreeBohr_to_Gpa +! write(u,"(I6,' ',13(2X,EN16.7))") i,i*sim%ts,sim%stat%et(i)*lo_Hartree_to_eV,sim%stat%ep(i)*lo_Hartree_to_eV,& +! sim%stat%ek(i)*lo_Hartree_to_eV,sim%stat%t(i),sim%stat%p(i)*lo_pressure_HartreeBohr_to_Gpa,sigma +! enddo +! close(u) +! +! write(*,*) '... wrote new stat' +! +! if ( allocated(sim%lattice) ) then +! u=open_file('out','outfile.lattice') +! do i=1,sim%nt +! do j=1,3 +! write(u,*) sim%lattice(:,j,i)*lo_Bohr_to_A +! enddo +! enddo +! close(u) +! write(*,*) '... wrote new lattice' +! endif + +end block init + +! check: block +! type(lo_symset) :: sym +! real(flyt), dimension(:,:), allocatable :: r0,r1,r2 +! real(flyt), dimension(:), allocatable :: lv0,lv1,x0,x1 +! real(flyt), dimension(3) :: v0,v1 +! real(flyt) :: f0,f1,f2 +! integer, dimension(:), allocatable :: di +! integer :: i,j,k,l,o,ii,jj +! integer :: ix +! +! ! call sym%generate(basis=uc%latticevectors,timereversal=.true.,verbosity=1) +! +! +! lo_allocate(r0(3,uc%na)) +! lo_allocate(di(uc%na)) +! r0=lo_huge +! di=1 +! do i=1,uc%na +! if ( di(i) .eq. 0 ) cycle +! ! round the coordinates of this atom +! do j=1,3 +! v0(j)=closest_fraction(uc%r(j,i)) +! enddo +! write(*,*) i,v0 +! do o=1,uc%sym%n +! v1=lo_operate_on_vector(uc%sym%op(o),v0,fractional=.true.) +! r0(:,uc%sym%op(o)%fmap(i))=v1 +! di(uc%sym%op(o)%fmap(i))=0 +! write(*,*) i,o,v1 +! enddo +! enddo +! +! write(*,*) +! do i=1,uc%na +! write(*,*) r0(:,i) +! enddo +! +! +! +! stop +! lo_allocate(r0(3,uc%na)) +! lo_allocate(r1(3,uc%na)) +! lo_allocate(r2(3,uc%na)) +! lo_allocate(lv0(3*uc%na)) +! lo_allocate(x0(ic%nx_internal)) +! lo_allocate(x1(ic%nx_internal)) +! ! Which degree of freedom? +! ix=1 +! ! Get the displacement pattern for this degree of freedom +! do i=1,uc%na +! lv0( (i-1)*3+1:i*3 )=uc%r(:,i) +! enddo +! x0=lo_clean_fractional_coordinates(matmul(transpose(ic%coeffM_internal),lv0) ) +! x0=matmul(transpose(ic%coeffM_internal),lv0) +! do i=1,ic%nx_internal +! x1(i)=closest_fraction(x0(i)) +! write(*,*) i,x0(i),x1(i) +! enddo +! +! lv0=lo_clean_fractional_coordinates(matmul(ic%coeffM_internal,x1)) +! lv0=matmul(ic%coeffM_internal,x0) +! +! do i=1,uc%na +! r1(:,i)=lv0( (i-1)*3+1:i*3 ) +! write(*,*) i +! write(*,*) uc%r(:,i) +! write(*,*) r1(:,i) +! write(*,*) uc%r(:,i)+r1(:,i) +! write(*,*) uc%r(:,i)-r1(:,i) +! enddo +! +! !write(*,*) shape(ic%coeffM_ +! +! +! write(*,*) matmul(transpose(ic%coeffM_internal),lv0) +! +! ! See if there is a nice one that fits! +! f0=lo_huge +! do ii=0,6 +! do jj=1,6 +! +! enddo +! enddo +! +! end block check + +contains + +function closest_fraction(x) result(y) + real(flyt) :: x, y + + integer, parameter :: nmax = 4 + real(flyt) :: f0, f1 + integer :: i, j + f0 = lo_huge + + do i = 0, nmax + do j = 1, nmax + f1 = real(i, flyt)/real(j, flyt) + if (abs(x - f1) .lt. f0) then + f0 = abs(x - f1) + y = f1 + end if + end do + end do +! y=x +end function + +end program + diff --git a/src/average_structure/options.f90 b/src/average_structure/options.f90 new file mode 100644 index 00000000..294feb1a --- /dev/null +++ b/src/average_structure/options.f90 @@ -0,0 +1,69 @@ +#include "precompilerdefinitions" +module options +use konstanter, only: flyt, lo_status, lo_author, lo_version, lo_licence, lo_hugeint +use flap, only: command_line_interface +implicit none + +private +public :: lo_opts + +type lo_opts + !> how much to talk + integer :: verbosity = -lo_hugeint + !> use every N timesteps + integer :: stride = -lo_hugeint + !> is it a constant pressure simulation? + logical :: npt = .false. +contains + procedure :: parse +end type + +contains + +subroutine parse(opts) + !> the options + class(lo_opts), intent(out) :: opts + !> the helper parser + type(command_line_interface) :: cli + ! + integer :: i + logical :: dumlog + real(flyt), dimension(3) :: dumflytv + + call cli%init(progname='average_structure', & + authors=lo_author, & + version=lo_version, & + license=lo_licence, & + help='Usage: ', & + description='Utility to fit equations of state to volume-energy data.', & + examples=["eosfit"], & + epilog=new_line('a')//"...") + + cli_manpage + cli_verbose + call cli%add(switch='--stride', switch_ab='-s', & + help='Number of dimenions for the fit', & + required=.false., act='store', def='1', error=lo_status) + if (lo_status .ne. 0) stop + + ! actually parse it + call cli%parse(error=lo_status) + if (lo_status .ne. 0) stop + + ! generate manpage? + call cli%get(switch='--manpage', val=dumlog) + if (dumlog) then + call cli%save_man_page(trim(cli%progname)//'.1') + call cli%save_usage_to_markdown(trim(cli%progname)//'.md') + write (*, *) 'Wrote manpage for "'//trim(cli%progname)//'"' + stop + end if + opts%verbosity = 0 + call cli%get(switch='--verbose', val=dumlog) + if (dumlog) opts%verbosity = 2 + + ! Parse the rest + call cli%get(switch='--stride', val=opts%stride, error=lo_status) +end subroutine + +end module diff --git a/src/average_structure/type_irredcell.f90 b/src/average_structure/type_irredcell.f90 new file mode 100644 index 00000000..74edaaf1 --- /dev/null +++ b/src/average_structure/type_irredcell.f90 @@ -0,0 +1,1031 @@ +#include "precompilerdefinitions" +!> Get the symmetry operations for a structure in a neat way. +module type_irredcell +! My stuff +use konstanter, only: flyt, lo_tol, lo_sqtol, lo_pi, lo_huge, lo_hugeint, lo_exitcode_symmetry, lo_status, & + lo_pressure_HartreeBohr_to_GPa, lo_force_HartreeBohr_to_eVA, & + lo_pressure_GPa_to_HartreeBohr, lo_pressure_eVA_to_HartreeBohr, & + lo_pressure_HartreeBohr_to_eVA, lo_Hartree_to_eV, lo_eV_to_Hartree, & + lo_bohr_to_A, lo_time_au_to_fs +use dump_data +use geometryfunctions +use gottochblandat !, only : lo_identitymatrix +use type_crystalstructure, only: lo_crystalstructure +use type_blas_lapack_wrappers, only: lo_dgesvd, lo_dgels, lo_dgglse, lo_dgelss, lo_dsyev +use type_mdsim, only: lo_mdsim +use type_symmetryoperation +implicit none + +private +public :: lo_irredcell +!public :: irreducible_trajectory +public :: avg_structure +public :: avg_structure_nvt +public :: remove_rotations + +!> symmetry-reduced representation of the unit cell +type lo_irredcell + !> number of degrees of freedom of the lattice + integer :: nx_lattice = -lo_hugeint + !> number of degrees of freedom for the forces + integer :: nx_internal = -lo_hugeint + !> coefficient matrix for the lattice + real(flyt), dimension(:, :), allocatable :: coeffM_lattice + !> coefficient matrix for the forces + real(flyt), dimension(:, :), allocatable :: coeffM_internal + !> reference lattice + real(flyt), dimension(3, 3) :: reference_lattice = -lo_huge + !> reference positions + real(flyt), dimension(:, :), allocatable :: reference_positions + !> irrep for the lattice + real(flyt), dimension(:), allocatable :: x_lattice + !> irrep for the positions + real(flyt), dimension(:), allocatable :: x_internal +contains + !> create everything + procedure :: generate + !> transform stress + !procedure :: strain_to_irreducible + + !> Store a configuration in the trajectory + !procedure :: add_step + !> Guess new structure + !procedure :: predict_new_structure + !> symmetrize stresses and forces + !procedure :: symmetrize_stress_and_forces +end type +contains + +! More stuff +#include "type_irredcell_aux.f90" + +!> Get the average of the internal structure +subroutine avg_structure_nvt(uc, ss, ic, sim, nuc, nss, ur) + !> unitcell + type(lo_crystalstructure), intent(in) :: uc + !> supercell + type(lo_crystalstructure), intent(in) :: ss + !> irreducible cell + type(lo_irredcell), intent(inout) :: ic + !> simulation + type(lo_mdsim), intent(inout) :: sim + !> averaged unitcell + type(lo_crystalstructure), intent(out) :: nuc + !> averaged supercell + type(lo_crystalstructure), intent(out) :: nss + !> Reference unitcell + type(lo_crystalstructure), intent(in), optional :: ur + + real(flyt), dimension(3, uc%na) :: average_positions + real(flyt), dimension(3, 3) :: average_lattice + + ! Align the simulation with the reference positions + init: block + real(flyt), dimension(:, :), allocatable :: dr + real(flyt), dimension(3*uc%na) :: flatr + real(flyt), dimension(3, 3) :: m0 + real(flyt), dimension(3) :: v0, v1 + real(flyt), dimension(ic%nx_internal) :: xavg + real(flyt) :: t0 + integer :: i + write (*, *) '' + write (*, *) 'Averaging structure' + + lo_allocate(dr(ic%nx_internal + 1, sim%nt)) + dr = 0.0_flyt + + t0 = walltime() + ! Align the simulation with the reference positions? Good idea? Not sure + xavg = 0.0_flyt + do i = 1, sim%nt + call get_internal_irrep_from_supercell(ic, uc, ss, sim%r(:, :, i)) + xavg = xavg + ic%x_internal + dr(1, i) = (i - 1)*sim%timestep*lo_time_au_to_fs + dr(2:1 + ic%nx_internal, i) = ic%x_internal + end do + xavg = xavg/sim%nt + + flatr = matmul(ic%coeffM_internal, xavg) + do i = 1, uc%na + average_positions(:, i) = ic%reference_positions(:, i) + flatr((i - 1)*3 + 1:i*3) + end do + average_lattice = uc%latticevectors + + call lo_dump_gnuplot_2d_real(dr, 'outfile.internal') + + write (*, *) 'got pos', walltime() - t0 + end block init + + ! And return clean proper structures + genclean: block + real(flyt), dimension(:, :), allocatable :: dr + real(flyt), dimension(3, 3) :: sslv, isslv + real(flyt), dimension(3) :: v0 + integer, dimension(3) :: ssd + integer :: i, j, l + + call nuc%generate(average_lattice, average_positions, uc%atomic_number, enhet=2) + call nuc%classify('bravais') + + ! Then the supercell. Have to be careful so that it is identical to the generating + ! supercell, so I will build it a little differently. + sslv = matmul(nuc%latticevectors, ss%info%supercellmatrix) + isslv = lo_invert3x3matrix(sslv) + + lo_allocate(dr(3, ss%na)) + do i = 1, ss%na + v0 = matmul(nuc%latticevectors, ss%info%cellindex(:, i) - 1) + nuc%rcart(:, ss%info%index_in_unitcell(i)) + dr(:, i) = matmul(isslv, v0) + end do + call nss%generate(sslv, dr, ss%atomic_number, enhet=2) + call nss%classify('supercell', nuc) + end block genclean + +end subroutine + +! !> Calculate the average structure from a simulation, in a neat symmetrized way. +! subroutine avg_internal_structure(uc,ss,ic,sim,nuc,nss) +! !> unitcell +! type(lo_crystalstructure), intent(in) :: uc +! !> supercell +! type(lo_crystalstructure), intent(in) :: ss +! !> irreducible cell +! type(lo_irredcell), intent(inout) :: ic +! !> simulation +! type(lo_mdsim), intent(inout) :: sim +! !> averaged unitcell +! type(lo_crystalstructure), intent(out) :: nuc +! !> averaged supercell +! type(lo_crystalstructure), intent(out) :: nss +! +! real(flyt), dimension(3,uc%na) :: average_positions +! real(flyt), dimension(3,3) :: average_lattice +! +! ! Get the average lattice +! if ( allocated(sim%lattice) ) then +! latt: block +! real(flyt), dimension(:,:), allocatable :: coeffM +! real(flyt), dimension(:), allocatable :: etaM +! real(flyt), dimension(3,3) :: I3,trf1,etaguess,m0,lattice +! real(flyt) :: f0,f1 +! integer :: t,i +! ! Some transformations +! I3=0.0_flyt +! do i=1,3 +! I3(i,i)=1.0_flyt +! enddo +! trf1=lo_invert3x3matrix( real(ss%info%supercellmatrix,flyt) ) +! ! First guess for eta? +! lo_allocate(coeffM(sim%nt*9,ic%nx_lattice)) +! lo_allocate(etaM(sim%nt*9)) +! coeffM=0.0_flyt +! etaM=0.0_flyt +! +! f0=0.0_flyt +! do t=1,sim%nt +! ! Current lattice +! lattice=matmul(sim%lattice(:,:,t),trf1) +! f0=f0+lo_determ(lattice) +! ! Convert to a strain and store for solution +! m0=matmul(lattice,uc%inv_latticevectors)-I3 +! coeffm( (t-1)*9+1:t*9,:)=ic%coeffM_lattice +! etam( (t-1)*9+1:t*9)=lo_flattentensor(m0) +! enddo +! f0=f0/sim%nt +! ! Solve +! call lo_linear_least_squares(coeffm,etam,ic%x_lattice) +! ! Multiply out to get the average +! m0=lo_unflatten_2tensor(matmul(ic%coeffM_lattice,ic%x_lattice))+I3 +! average_lattice=matmul(m0,ic%reference_lattice) +! f1=(f0/lo_determ(average_lattice))**(1.0_flyt/3.0_flyt) +! average_lattice=average_lattice*f1 +! +! write(*,*) 'latvec',lo_determ(average_lattice),f1,f0 +! do i=1,3 +! write(*,*) average_lattice(:,i)*lo_bohr_to_A +! enddo +! end block latt +! endif +! +! ! And the positions. Maybe a little trickier, maybe not. +! posi: block +! integer, parameter :: maxnavg=10 +! integer, parameter :: stride=2 +! real(flyt), dimension(:,:), allocatable :: ucref,ssref,ucpos,sspos,r0,r1 +! real(flyt), dimension(:,:), allocatable :: temp_avgpos_uc,temp_refpos_uc,refpos +! +! +! real(flyt), dimension(:,:), allocatable :: irr_internal,dr,du +! real(flyt), dimension(:,:), allocatable :: wA +! real(flyt), dimension(:), allocatable :: WB +! real(flyt), dimension(3,3) :: trf1,trf2,lattice,invlattice +! real(flyt), dimension(3) :: v0,v1,v2,refshift +! real(flyt) :: f0,f1,f2,f3 +! integer, dimension(:), allocatable :: dctr +! integer :: nf,nt,tt,t,i,j,k,ii,jj +! integer :: iter +! +! ! Space for strange averages and things. +! lo_allocate(sspos(3,ss%na)) +! lo_allocate(ssref(3,ss%na)) +! lo_allocate(ucref(3,uc%na)) +! lo_allocate(ucpos(3,uc%na)) +! lo_allocate(r0(3,uc%na)) +! lo_allocate(r1(3,uc%na)) +! ucref=uc%r +! ssref=ss%r +! +! ! Get unwrapped positions +! do t=1,sim%nt +! do i=1,sim%na +! v0=lo_clean_fractional_coordinates(sim%r(:,i,t)-ssref(:,i)+0.5_flyt)-0.5_flyt +! sim%r_npbc(:,i,t)=ssref(:,i)+v0 +! enddo +! enddo +! +! ! Raw average positions for the supercell +! sspos=0.0_flyt +! do t=1,sim%nt +! do i=1,sim%na +! sspos(:,i)=sspos(:,i)+sim%r_npbc(:,i,t) +! enddo +! enddo +! sspos=sspos/sim%nt +! +! ! Downfolded average positions +! f0=uc%volume/ss%volume +! trf1=matmul(uc%inv_latticevectors,ss%latticevectors) +! ucpos=0.0_flyt +! do i=1,ss%na +! j=ss%info%index_in_unitcell(i) +! v0=ss%info%cellindex(:,i)-1 +! v0=matmul(trf1,sspos(:,i))-v0 +! ucpos(:,j)=ucpos(:,j)+v0*f0 +! enddo +! +! ! Symmetrize these? +! r0=ucpos +! r1=0.0_flyt +! do iter=1,2 +! do ii=1,uc%sym%n +! do i=1,uc%na +! j=uc%sym%op(ii)%fmap(i) +! r1(:,j)=lo_operate_on_vector(uc%sym%op(ii),r0(:,i),fractional=.true.) +! enddo +! r1=lo_clean_fractional_coordinates(r1-r0+0.5_flyt,1E-13_flyt)-0.5_flyt +! r0=lo_clean_fractional_coordinates(r0+0.5_flyt*r1,1E-13_flyt) +! enddo +! r0=lo_chop(r0,1E-13_flyt) +! enddo +! +! average_positions=r0 +! ! do i=1,uc%na +! ! average_positions(:,i)=r0(:,i)-r0(:,1) +! ! enddo +! end block posi +! +! ! And return clean proper structures +! genclean: block +! real(flyt), dimension(:,:), allocatable :: dr +! real(flyt), dimension(3,3) :: sslv,isslv +! real(flyt), dimension(3) :: v0 +! integer, dimension(3) :: ssd +! integer :: i,j,l +! +! call nuc%generate(average_lattice,average_positions,uc%atomic_number,enhet=2) +! call nuc%classify('bravais') +! +! ! Then the supercell. Have to be careful so that it is identical to the generating +! ! supercell, so I will build it a little differently. +! sslv=matmul(nuc%latticevectors,ss%info%supercellmatrix) +! isslv=lo_invert3x3matrix(sslv) +! +! lo_allocate(dr(3,ss%na)) +! do i=1,ss%na +! v0=matmul(nuc%latticevectors,ss%info%cellindex(:,i)-1)+nuc%rcart(:,ss%info%index_in_unitcell(i)) +! dr(:,i)=matmul(isslv,v0) +! enddo +! call nss%generate(sslv,dr,ss%atomic_number,enhet=2) +! call nss%classify('supercell',nuc) +! ! +! end block genclean +! +! end subroutine + +!> Remove arbitrary precession of the lattice +subroutine remove_rotations(uc, ss, ic, sim) + !> unitcell + type(lo_crystalstructure), intent(in) :: uc + !> supercell + type(lo_crystalstructure), intent(in) :: ss + !> irreducible cell + type(lo_irredcell), intent(inout) :: ic + !> simulation + type(lo_mdsim), intent(inout) :: sim + + real(flyt), dimension(3, 3) :: rotation, etaguess, lattice, strain + real(flyt), dimension(3, 3) :: m0, invlatt, trf1, trf2, I3 + real(flyt), dimension(3) :: v0, v1, com, v2 + real(flyt) :: invmass, invna + integer :: i, t + + ! Some transformations + I3 = 0.0_flyt + do i = 1, 3 + I3(i, i) = 1.0_flyt + end do + trf1 = lo_invert3x3matrix(real(ss%info%supercellmatrix, flyt)) + trf2 = real(ss%info%supercellmatrix, flyt) + ! First guess for eta? + m0 = matmul(sim%lattice(:, :, 1), trf1) + m0 = matmul(m0, invlatt) + etaguess = m0 + + ! First, temporarily convert everything to fractional coordinates. + do t = 1, sim%nt + lattice = matmul(sim%lattice(:, :, t), trf1) + invlatt = lo_invert3x3matrix(lattice) + ! Forces and velocities + do i = 1, sim%na + sim%f(:, i, t) = matmul(invlatt, sim%f(:, i, t)) + sim%v(:, i, t) = matmul(invlatt, sim%v(:, i, t)) + end do + ! Set things to zero that have to be fixed later + sim%u(:, :, t) = 0.0_flyt + ! And the stress tensor. Hmm. I think it works like this. + ! Actually not sure. + m0 = sim%stat%stress(:, :, t) + m0 = matmul(matmul(invlatt, m0), transpose(invlatt)) + sim%stat%stress(:, :, t) = m0 + end do + + ! First guess for eta? + m0 = matmul(sim%lattice(:, :, 1), trf1) + m0 = matmul(m0, invlatt) + etaguess = m0 + ! Now find the rotaiton for the lattice, and rotate it + do t = 1, sim%nt + ! Get the transform for the lattice + lattice = matmul(sim%lattice(:, :, t), trf1) + call lo_find_rotation_that_makes_strain_diagonal(ic%reference_lattice, lattice, rotation, strain, etaguess) + etaguess = strain + ! Then transform it and store + m0 = matmul(rotation, lattice) + m0 = matmul(m0, trf2) + sim%lattice(:, :, t) = m0 + end do + + ! Then convert everything back to Cartesian coordinates + do t = 1, sim%nt + lattice = matmul(sim%lattice(:, :, t), trf1) + ! And the stress tensor. Hmm. I think it works like this. + m0 = sim%stat%stress(:, :, t) + m0 = matmul(matmul(lattice, m0), transpose(lattice)) + sim%stat%stress(:, :, t) = m0 + ! Forces, velocities and displacements. Not sure what makes more + ! sense, to use some average or the instantaneous. Will use the + ! instantaneous, I concluded for no good reason. + do i = 1, sim%na + sim%f(:, i, t) = matmul(lattice, sim%f(:, i, t)) + sim%v(:, i, t) = matmul(lattice, sim%v(:, i, t)) + v0 = lo_clean_fractional_coordinates(sim%r(:, i, t) - ss%r(:, i) + 0.5_flyt) - 0.5_flyt + sim%r_npbc(:, i, t) = sim%r(:, i, t) + v0 + sim%u(:, i, t) = matmul(lattice, v0) + end do + end do + + invmass = 1.0_flyt/sum(ss%mass) + invna = 1.0_flyt/real(ss%na, flyt) + + ! In the new coordinates, also remove any drift. + do t = 1, sim%nt + v0 = 0.0_flyt + v1 = 0.0_flyt + do i = 1, sim%na + v2 = lo_clean_fractional_coordinates(sim%r(:, i, t) - ss%r(:, i) + 0.5_flyt) - 0.5_flyt + v0 = v0 + v2 !*ss%mass(i) + v1 = v1 + sim%f(:, i, t) + end do + + v0 = v0*invna !mass + v1 = v1*invna + do i = 1, sim%na + sim%r(:, i, t) = sim%r(:, i, t) - v0 + sim%f(:, i, t) = sim%f(:, i, t) - v1 + v2 = lo_clean_fractional_coordinates(sim%r(:, i, t) - ss%r(:, i) + 0.5_flyt) - 0.5_flyt + sim%r_npbc(:, i, t) = sim%r(:, i, t) + v2 + sim%u(:, i, t) = matmul(sim%lattice(:, :, t), v2) + end do + end do +end subroutine + +!> Calculate the average structure from a simulation, in a neat symmetrized way. +subroutine avg_structure(uc, ss, ic, sim, nuc, nss) + !> unitcell + type(lo_crystalstructure), intent(in) :: uc + !> supercell + type(lo_crystalstructure), intent(in) :: ss + !> irreducible cell + type(lo_irredcell), intent(inout) :: ic + !> simulation + type(lo_mdsim), intent(inout) :: sim + !> averaged unitcell + type(lo_crystalstructure), intent(out) :: nuc + !> averaged supercell + type(lo_crystalstructure), intent(out) :: nss + + real(flyt), dimension(3, uc%na) :: average_positions + real(flyt), dimension(3, 3) :: average_lattice + + ! Get the average lattice + if (allocated(sim%lattice)) then + latt: block + real(flyt), dimension(:, :), allocatable :: coeffM + real(flyt), dimension(:), allocatable :: etaM + real(flyt), dimension(3, 3) :: I3, trf1, etaguess, m0, lattice + real(flyt) :: f0, f1 + integer :: t, i + ! Some transformations + I3 = 0.0_flyt + do i = 1, 3 + I3(i, i) = 1.0_flyt + end do + trf1 = lo_invert3x3matrix(real(ss%info%supercellmatrix, flyt)) + ! First guess for eta? + lo_allocate(coeffM(sim%nt*9, ic%nx_lattice)) + lo_allocate(etaM(sim%nt*9)) + coeffM = 0.0_flyt + etaM = 0.0_flyt + + f0 = 0.0_flyt + do t = 1, sim%nt + ! Current lattice + lattice = matmul(sim%lattice(:, :, t), trf1) + f0 = f0 + lo_determ(lattice) + ! Convert to a strain and store for solution + m0 = matmul(lattice, uc%inv_latticevectors) - I3 + coeffm((t - 1)*9 + 1:t*9, :) = ic%coeffM_lattice + etam((t - 1)*9 + 1:t*9) = lo_flattentensor(m0) + end do + f0 = f0/sim%nt + ! Solve + call lo_linear_least_squares(coeffm, etam, ic%x_lattice) + ! Multiply out to get the average + m0 = lo_unflatten_2tensor(matmul(ic%coeffM_lattice, ic%x_lattice)) + I3 + average_lattice = matmul(m0, ic%reference_lattice) + f1 = (f0/lo_determ(average_lattice))**(1.0_flyt/3.0_flyt) + average_lattice = average_lattice*f1 + + write (*, *) 'latvec', lo_determ(average_lattice), f1, f0 + do i = 1, 3 + write (*, *) average_lattice(:, i)*lo_bohr_to_A + end do + end block latt + else + average_lattice = uc%latticevectors + end if + + ! And the positions. Maybe a little trickier, maybe not. + posi: block + integer, parameter :: maxnavg = 10 + integer, parameter :: stride = 2 + real(flyt), dimension(:, :), allocatable :: ucref, ssref, ucpos, sspos, r0, r1 + real(flyt), dimension(:, :), allocatable :: temp_avgpos_uc, temp_refpos_uc, refpos + + real(flyt), dimension(:, :), allocatable :: irr_internal, dr, du + real(flyt), dimension(:, :), allocatable :: wA + real(flyt), dimension(:), allocatable :: WB + real(flyt), dimension(3, 3) :: trf1, trf2, lattice, invlattice + real(flyt), dimension(3) :: v0, v1, v2, refshift + real(flyt) :: f0, f1, f2, f3 + integer, dimension(:), allocatable :: dctr + integer :: nf, nt, tt, t, i, j, k, ii, jj + integer :: iter + + ! Space for strange averages and things. + lo_allocate(sspos(3, ss%na)) + lo_allocate(ssref(3, ss%na)) + lo_allocate(ucref(3, uc%na)) + lo_allocate(ucpos(3, uc%na)) + lo_allocate(r0(3, uc%na)) + lo_allocate(r1(3, uc%na)) + ucref = uc%r + ssref = ss%r + + ! Get unwrapped positions + do t = 1, sim%nt + do i = 1, sim%na + v0 = lo_clean_fractional_coordinates(sim%r(:, i, t) - ssref(:, i) + 0.5_flyt) - 0.5_flyt + sim%r_npbc(:, i, t) = ssref(:, i) + v0 + end do + end do + write (*, *) sum(abs(sim%r_npbc)) + + ! Raw average positions for the supercell + sspos = 0.0_flyt + do t = 1, sim%nt + do i = 1, sim%na + sspos(:, i) = sspos(:, i) + sim%r_npbc(:, i, t) + end do + end do + sspos = sspos/sim%nt + + ! Downfolded average positions + f0 = uc%volume/ss%volume + trf1 = matmul(uc%inv_latticevectors, ss%latticevectors) + ucpos = 0.0_flyt + do i = 1, ss%na + j = ss%info%index_in_unitcell(i) + v0 = ss%info%cellindex(:, i) - 1 + v0 = matmul(trf1, sspos(:, i)) - v0 + ucpos(:, j) = ucpos(:, j) + v0*f0 + end do + + ! Symmetrize these? + r0 = ucpos + do iter = 1, 2 + do ii = 1, uc%sym%n + r1 = 0.0_flyt + do i = 1, uc%na + j = uc%sym%op(ii)%fmap(i) + r1(:, j) = lo_operate_on_vector(uc%sym%op(ii), r0(:, i), fractional=.true.) + end do + r1 = lo_clean_fractional_coordinates(r1 - r0 + 0.5_flyt, 1E-13_flyt) - 0.5_flyt + r0 = lo_clean_fractional_coordinates(r0 + 0.5_flyt*r1, 1E-13_flyt) + end do + r0 = lo_chop(r0, 1E-13_flyt) + end do + + average_positions = r0 +! do i=1,uc%na +! average_positions(:,i)=r0(:,i)-r0(:,1) +! enddo + end block posi + + ! And return clean proper structures + genclean: block + real(flyt), dimension(:, :), allocatable :: dr + real(flyt), dimension(3, 3) :: sslv, isslv + real(flyt), dimension(3) :: v0 + integer, dimension(3) :: ssd + integer :: i, j, l + + call nuc%generate(average_lattice, average_positions, uc%atomic_number, enhet=2) + call nuc%classify('bravais') + + ! Then the supercell. Have to be careful so that it is identical to the generating + ! supercell, so I will build it a little differently. + sslv = matmul(nuc%latticevectors, ss%info%supercellmatrix) + isslv = lo_invert3x3matrix(sslv) + + lo_allocate(dr(3, ss%na)) + do i = 1, ss%na + v0 = matmul(nuc%latticevectors, ss%info%cellindex(:, i) - 1) + nuc%rcart(:, ss%info%index_in_unitcell(i)) + dr(:, i) = matmul(isslv, v0) + end do + call nss%generate(sslv, dr, ss%atomic_number, enhet=2) + call nss%classify('supercell', nuc) + ! + end block genclean + +end subroutine + +!!> Get the irreducible trajectory +!subroutine irreducible_trajectory(uc,ss,ic,sim) +! !> unitcell +! type(lo_crystalstructure), intent(in) :: uc +! !> supercell +! type(lo_crystalstructure), intent(in) :: ss +! !> irreducible cell +! type(lo_irredcell), intent(inout) :: ic +! !> simulation +! type(lo_mdsim), intent(inout) :: sim +! +! real(flyt), dimension(3,3) :: avg_lattice +! +! +! +! +! ! First the lattice +! latt: block +! real(flyt), dimension(:,:), allocatable :: irr_lattice +! real(flyt), dimension(9) :: fv0 +! real(flyt), dimension(3,3) :: m0,trf1,I3 +! integer :: i +!real(flyt) :: f0,f1,f2 +! +! I3=0.0_flyt +! do i=1,3 +! I3(i,i)=1.0_flyt +! enddo +! +! ! Transformation matrix to unitcell +! trf1=lo_invert3x3matrix( real(ss%info%supercellmatrix,flyt) ) +! lo_allocate(irr_lattice(ic%nx_lattice,sim%nt)) +! irr_lattice=0.0_flyt +! +! f0=0.0_flyt +! do i=1,sim%nt +! ! Get the unitcell strain. Always useful. First convert to unitcell lattice vectors +! m0=matmul(trf1,sim%lattice(:,:,i)) +! call ic%strain_to_irreducible(m0) +! irr_lattice(:,i)=ic%x_lattice +! f0=f0+lo_determ(sim%lattice(:,:,i))/sim%na/sim%nt +! enddo +! write(*,*) 'avg vol/A',f0 +! do i=1,ic%nx_lattice +! ic%x_lattice(i)=lo_mean(irr_lattice(i,:)) +! enddo +! fv0=matmul(ic%coeffM_lattice,ic%x_lattice) +! avg_lattice=matmul(ic%reference_lattice,I3+lo_unflatten_2tensor(fv0)) +! +!f0=lo_determ(avg_lattice)/uc%na +!write(*,*) 'avg vol/A',f0 +!! +!!fv0=matmul(ic%coeffM_lattice,[1.0_flyt,0.5_flyt]) +!!avg_lattice=lo_unflatten_2tensor(fv0) +!! +!m0=lo_unflatten_2tensor(fv0) +! +!do i=1,3 +! write(*,*) m0(:,i) !avg_lattice(:,i)*lo_bohr_to_A +!enddo +!do i=1,3 +! write(*,*) ic%reference_lattice(:,i) +!enddo +!do i=1,3 +! write(*,*) avg_lattice(:,i)/norm2(avg_lattice(:,i)) +!enddo +!stop +! +!! +! +! end block latt +! +! ! Think about center of mass, drift and such. +! centerofmass: block +! real(flyt), dimension(3) :: v0,v1,v2 +! integer :: i,j,k,l +! +! ! First align the centers of mass with each other. +! ! do i=1,sim%nt +! ! v0=0.0_flyt +! ! do j=1,sim%na +! ! v1=sim%r_npbc(:,j,i)-ss%r(:,j) +! ! v0=v0+v1*ss%mass(j) +! ! enddo +! ! v0=v0/sum(ss%mass) +! ! do j=1,sim%na +! ! sim%r(:,j,i)=sim%r(:,j,i)-v0 +! ! sim%r_npbc(:,j,i)=sim%r_npbc(:,j,i)-v0 +! ! enddo +! ! enddo +! ! ! Then make sure there is no rigid shift of everything. +! ! v0=0.0_flyt +! ! do i=1,sim%nt +! ! do j=1,sim%na +! ! v1=lo_clean_fractional_coordinates(sim%r_npbc(:,j,i)-ss%r(:,j)+0.5_flyt)-0.5_flyt +! ! v0=v0+v1 +! ! enddo +! ! enddo +! ! v0=v0/sim%nt/sim%na +! ! do i=1,sim%nt +! ! do j=1,sim%na +! ! sim%r(:,j,i)=sim%r(:,j,i)-v0 +! ! sim%r_npbc(:,j,i)=sim%r_npbc(:,j,i)-v0 +! ! enddo +! ! enddo +! end block centerofmass +! +! ! Then the internal positions. Hmm. +! posi: block +! integer, parameter :: maxnavg=10 +! real(flyt), dimension(:,:), allocatable :: irr_internal +! real(flyt), dimension(:,:), allocatable :: big_coeff,big_ref,wA,WB +! real(flyt), dimension(:), allocatable :: big_pos +! real(flyt), dimension(3,3) :: trf1 +! real(flyt), dimension(3) :: v0,v1 +! real(flyt), dimension(:), allocatable :: invctr +! integer :: i,j,k,l,ii +! integer :: navg,nf +! +! ! How many bins for the averaging? +! navg=min(sim%nt,maxnavg) +! nf=sim%na*3 +! +! ! First get the reference positions +! lo_allocate(big_ref(3,ss%na)) +! lo_allocate(big_pos(3*ss%na)) +! lo_allocate(big_coeff(3*ss%na,ic%nx_internal)) +! lo_allocate(irr_internal(ic%nx_internal,sim%nt)) +! irr_internal=0.0_flyt +! big_ref=0.0_flyt +! +! ! Get coefficient matrix and reference positions +! big_coeff=0.0_flyt +! do i=1,ss%na +! j=ss%info%index_in_unitcell(i) +! v0=ic%reference_positions(:,j) +! v0=v0+(ss%info%cellindex(:,i)-1) +! big_ref( :,i )=v0 +! big_coeff( (i-1)*3+1:i*3,: )=ic%coeffM_internal( (j-1)*3+1:j*3,:) +! enddo +! ! Expand the big matrix +! lo_allocate(wA(nf*navg,ic%nx_internal)) +! lo_allocate(wB(nf*navg,1)) +! wA=0.0_flyt +! wB=0.0_flyt +! do i=1,navg +! wA( (i-1)*nf+1:i*nf,: )=big_coeff +! enddo +! +! ! Get the averaging factor +! lo_allocate(invctr(navg)) +! invctr=0.0_flyt +! do i=1,sim%nt +! j=mod(i,navg)+1 +! invctr(j)=invctr(j)+1.0_flyt +! enddo +! invctr=1.0_flyt/invctr +! +! trf1=real(ss%info%supercellmatrix,flyt) +! sim%r_npbc(:,:,1)=ss%r +! do i=1,sim%nt +! +! do j=1,sim%na +! v0=matmul(trf1,sim%r_npbc(:,j,i))-big_ref(:,j) +! big_pos( (j-1)*3+1:j*3 )=v0 +! enddo +! !big_pos=lo_clean_fractional_coordinates(big_pos+0.5_flyt)-0.5_flyt +! irr_internal(:,i)=matmul(transpose(big_coeff),big_pos) +!write(*,*) i,irr_internal(:,i) +! enddo +! +! do i=1,sim%na +! do j=1,3 +! big_ref(j,i)=lo_mean(sim%r_npbc(j,i,:)) +! enddo +! enddo +! big_ref=lo_clean_fractional_coordinates(big_ref) +!write(*,*) 'Average:' +! do i=1,sim%na +! write(*,*) i,lo_clean_fractional_coordinates(big_ref(:,i)-ss%r(:,i)+0.5_flyt)-0.5_flyt +! enddo +! +! end block posi +! +!end subroutine + +!> create the irreducible representation of the lattice +subroutine generate(ic, p, verbosity, reference_positions) + !> irreducible representation + class(lo_irredcell), intent(out) :: ic + !> normal crystal structure + type(lo_crystalstructure), intent(inout) :: p + !> talk? + integer, intent(in) :: verbosity + !> Not the normal reference positions? + real(flyt), dimension(:, :), intent(in), optional :: reference_positions + + ! decide on settings and such + init: block + integer :: i + if (verbosity .gt. 0) then + write (*, *) '' + write (*, *) 'BUILDING IRREDUCIBLE REPRESENTATION OF THE LATTICE' + end if + + ! We need all the symmetry things + if (p%info%havespacegroup .eqv. .false.) then + call p%classify('wedge', timereversal=.true.) + end if + + ! Store the reference positions + lo_allocate(ic%reference_positions(3, p%na)) + if (present(reference_positions)) then + ic%reference_positions = reference_positions + else + ic%reference_positions = p%r + end if + ic%reference_lattice = p%latticevectors + + end block init + + ! Start with the lattice vectors, or more precisely the symmetry of the strain tensor. + lattice: block + real(flyt), dimension(9, 9) :: invarM, I9, T, rotM + real(flyt), dimension(3, 3) :: m0 + real(flyt), dimension(9) :: dv + integer :: i, j + + ! 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_flyt + do i = 1, p%sym%n + rotM = lo_expandoperation_pair(p%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, ic%coeffM_lattice, ic%nx_lattice, 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 (ic%nx_lattice .lt. 1 .or. ic%nx_lattice .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 + lo_allocate(ic%x_lattice(ic%nx_lattice)) + ic%x_lattice = 0.0_flyt + + if (verbosity .gt. 0) then + write (*, *) '... found ', tochar(ic%nx_lattice), ' degrees of freedom for the lattice:' + do i = 1, ic%nx_lattice + write (*, *) 'x'//tochar(i)//':' + dv = 0.0_flyt + dv(i) = 1.0_flyt + m0 = lo_unflatten_2tensor(matmul(ic%coeffM_lattice, dv(1:ic%nx_lattice))) + do j = 1, 3 + write (*, "(3(2X,F16.9))") m0(:, j) + end do + end do + end if + + end block lattice + + ! Now do the same thing, but for the internal positions. And in this case it's the symmetry of the + ! forces we are interested in. + bigmatrix: block + real(flyt), dimension(:, :), allocatable :: invarM, rotM, IM, dr + real(flyt), dimension(:), allocatable :: dv + integer :: i, j, o, n, u + ! + n = p%na*3 + ! Some space for stuff + lo_allocate(rotM(n, n)) + lo_allocate(IM(n, n)) + lo_allocate(invarM(n, n)) + call lo_identitymatrix(IM) + + ! Build a matrix representation how all the internal positions transform under an + ! operation. I use large block matrices for this. + invarM = 0.0_flyt + do o = 1, p%sym%n + rotM = 0.0_flyt + do i = 1, p%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 = p%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(p%sym%op(o)%fm) + end do + invarM = invarM + rotM - IM + 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_flyt + do i = 1, p%na + do j = 1, 3 + rotM((i - 1)*3 + j, j) = 1.0_flyt + end do + end do + invarM = invarM + rotM + + ! In an ideal world I should probably add rotational invariance as well, but I don't think + ! it matters, and I also don't know how to do it. Also worth noting is that the symmetry is + ! done in fractional coordinates, and all shifts from starting positions are also in fractional + ! coordinates. Otherwise it gets real wonky when the volume changes. + + ! Anyway, now we can project out the nullspace of this matrix + call lo_nullspace_coefficient_matrix(invarM, ic%coeffM_internal, ic%nx_internal, tolerance=lo_sqtol) + ! and some space for a solution + if (ic%nx_internal .gt. 0) then + lo_allocate(ic%x_internal(ic%nx_internal)) + ic%x_internal = 0.0_flyt + end if + + if (verbosity .gt. 0) then + write (*, *) '... found ', tochar(ic%nx_internal), ' internal degrees of freedom' + lo_allocate(dv(p%na*3)) + lo_allocate(dr(3, p%na)) + do i = 1, ic%nx_internal + write (*, *) 'x'//tochar(i)//': (fractional, Cartesian)' + dv = 0.0_flyt + dv(i) = 1.0_flyt + dv = matmul(ic%coeffM_internal, dv(1:ic%nx_internal)) + do j = 1, p%na + dr(:, j) = dv((j - 1)*3 + 1:j*3) + write (*, "(6(2X,F15.9))") dr(:, j), matmul(p%latticevectors, dr(:, j)) + end do + end do + write (*, *) '' + write (*, *) 'Reference positions:' + do j = 1, p%na + dr(:, j) = dv((j - 1)*3 + 1:j*3) + write (*, "(3(2X,F15.9))") ic%reference_positions(:, j) + end do + end if + end block bigmatrix + +end subroutine + +! !> Turn the raw forces and stress tensor from AIMS into something clean and symmetric. Not used for anything at the moment. But could come handy on a rainy day. +! subroutine symmetrize_stress_and_forces(ic,p,orig_stress,orig_forces,clean_stress,clean_forces) +! !> irreducible representation +! class(lo_irredcell), intent(inout) :: ic +! !> crystal structure +! type(lo_crystalstructure), intent(in) :: p +! !> noisy stress tensor +! real(flyt), dimension(3,3), intent(in) :: orig_stress +! !> noisy forces +! real(flyt), dimension(:,:), intent(in) :: orig_forces +! !> nice stress tensor +! real(flyt), dimension(3,3), intent(out) :: clean_stress +! !> nice forces +! real(flyt), dimension(:,:), intent(out) :: clean_forces +! +! ! least square force the stress into the irreducible form +! fixlattice: block +! real(flyt), dimension(3,3) :: m0 +! real(flyt), dimension(9,ic%nx_lattice) :: wM +! real(flyt), dimension(9,1) :: wV +! real(flyt) :: f0 +! integer :: i +! +! ! Solve for the irreducible components +! wM=ic%coeffM_lattice +! wV(:,1)=lo_flatten_2tensor(orig_stress) +! call lo_dgels(wM,wV,info=lo_status) +! if ( lo_status .ne. 0 ) then +! write(*,*) 'Lapack failed for the lattice' +! stop +! endif +! ic%x_lattice=wV(1:ic%nx_lattice,1) +! ! Store the symetrized stress +! clean_stress=lo_chop( lo_unflatten_2tensor(matmul(ic%coeffM_lattice,ic%x_lattice)), 1E-12_flyt) +! end block fixlattice +! +! ! Same thing for the forces +! if ( ic%nx_internal .gt. 0 ) then +! fixforces: block +! real(flyt), dimension(:,:), allocatable :: mA,mC,mB +! real(flyt), dimension(:), allocatable :: vB,vD +! real(flyt), dimension(3) :: v0 +! integer :: i,j,l,n +! +! n=size(ic%coeffM_internal,1) +! if ( size(orig_forces,1)*size(orig_forces,2) .ne. n ) then +! write(*,*) 'Wrong number of forces' +! stop +! endif +! +! ! normal least squares +! lo_allocate(mA(n,ic%nx_internal)) +! lo_allocate(mB(n,1)) +! l=0 +! do i=1,size(orig_forces,2) +! do j=1,size(orig_forces,1) +! l=l+1 +! mB(l,1)=orig_forces(j,i) +! enddo +! enddo +! call lo_dgels(mA,mB,info=lo_status) +! if ( lo_status .ne. 0 ) then +! write(*,*) 'Lapack failed for forces' +! stop +! endif +! ic%x_internal=mB(1:ic%nx_internal,1) +! +! ! Store the symmetrized forces +! if ( .not.allocated(vB) ) allocate(vB(n)) +! vB=matmul(ic%coeffM_internal,ic%x_internal) +! do i=1,p%na +! v0=vb( (i-1)*3+1:i*3 ) +! v0=matmul( p%latticevectors, v0 ) +! clean_forces(:,i)=lo_chop( v0, 1E-12_flyt ) +! enddo +! end block fixforces +! else +! ! In case there are no internal degrees of freedom. +! clean_forces=0.0_flyt +! endif +! end subroutine + +end module diff --git a/src/average_structure/type_irredcell_aux.f90 b/src/average_structure/type_irredcell_aux.f90 new file mode 100644 index 00000000..413bd4a9 --- /dev/null +++ b/src/average_structure/type_irredcell_aux.f90 @@ -0,0 +1,430 @@ + +!> Get irreducible representation from something +subroutine get_internal_irrep_from_unitcell(ic, r) + !> irreducible representation + class(lo_irredcell), intent(inout) :: ic + !> unitcell positions in fractional coordinates + real(flyt), dimension(:, :), intent(in) :: r + + real(flyt), dimension(:), allocatable :: flatu + real(flyt), dimension(3) :: v0 + integer :: i + + lo_allocate(flatu(size(r))) + flatu = 0.0_flyt + do i = 1, size(r, 2) + v0 = lo_clean_fractional_coordinates(r(:, i) - ic%reference_positions(:, i) + 0.5_flyt) - 0.5_flyt + flatu((i - 1)*3 + 1:i*3) = v0 + end do + call lo_linear_least_squares(ic%coeffM_internal, flatu, ic%x_internal) +end subroutine + +!> Get irreducible representation +subroutine get_internal_irrep_from_supercell(ic, uc, ss, r) + !> irreducible representation + class(lo_irredcell), intent(inout) :: ic + !> unitcell + type(lo_crystalstructure), intent(in) :: uc + !> supercell + type(lo_crystalstructure), intent(in) :: ss + !> unitcell positions in fractional coordinates + real(flyt), dimension(:, :), intent(in) :: r + + real(flyt), dimension(:, :), allocatable :: coeffM + real(flyt), dimension(:), allocatable :: flatu + real(flyt), dimension(3, 3) :: m0 + real(flyt), dimension(3) :: v0 + integer :: i, j + + lo_allocate(flatu(size(r))) + lo_allocate(coeffM(size(r), ic%nx_internal)) + flatu = 0.0_flyt + coeffM = 0.0_flyt + m0 = matmul(uc%inv_latticevectors, ss%latticevectors) + + do i = 1, size(r, 2) + j = ss%info%index_in_unitcell(i) + v0 = matmul(m0, r(:, i)) - ss%info%cellindex(:, i) + 1 - ic%reference_positions(:, j) + v0 = lo_clean_fractional_coordinates(v0 + 0.5_flyt) - 0.5_flyt + flatu((i - 1)*3 + 1:i*3) = v0 + coeffM((i - 1)*3 + 1:i*3, :) = ic%coeffM_internal((j - 1)*3 + 1:j*3, :) + end do + call lo_linear_least_squares(coeffM, flatu, ic%x_internal) +end subroutine + +!!> Finds the rotation matrix such that lattice'=rotation*lattice represents a diagnoal distortion of ref_lattice. +!subroutine lo_find_rotation_that_makes_strain_diagonal( ref_lattice,lattice,rotation,strain,guess ) +! !> reference lattice +! real(flyt), dimension(3,3), intent(in) :: ref_lattice +! !> current lattice +! real(flyt), dimension(3,3), intent(in) :: lattice +! !> rotation matrix +! real(flyt), dimension(3,3), intent(out) :: rotation +! !> strain matrix +! real(flyt), dimension(3,3), intent(out) :: strain +! !> guess for the symmetric strain? +! real(flyt), dimension(3,3), intent(in), optional :: guess +! +! integer, parameter :: nouter=1000 +! integer, parameter :: ninner=1000 +! integer, parameter :: maxctr=20 +! real(flyt), parameter :: defstep=1E-5_flyt +! real(flyt), parameter :: errtol=1E-20_flyt +! +! real(flyt), dimension(6,6,6) :: coeffM +! real(flyt), dimension(6) :: bvec,eta0,eta1,grad +! real(flyt), dimension(3,3) :: m0 +! real(flyt) :: e0,e1,f0,step +! integer :: initer,outiter,ctr,i +! +! ! So, what do I solve here. The idea is that any transformation of a lattice takes +! ! a,b,c,alpha,beta,gamma to new values. I need six numbers to define that transformation +! ! uniquely. But you often come across these transformations as 3x3-matrices, so that they +! ! might contain an extra rotation. This routine finds the diagonal form of a general +! ! transformation between lattice vectors, since that can always be done. +! ! +! ! So, in practice, I write the transformed lattice vectors +! ! +! ! a'=eta*a +! ! b'=eta*b +! ! c'=eta*c +! ! +! ! And find the symmetric matrix eta that produces the correct values for +! ! +! ! a.a b.b c.c a.b a.c b.c +! ! +! ! when applied to the original lattice. I end up with a bunch of matrix +! ! equations on the form X^T A X = b, that I solve with a steepest descent. +! ! Seems to be a robust procedure. +! +! ! Get the coefficient matrices, A in XAX=b +! call cfvars(coeffM(:,:,1),coeffM(:,:,2),coeffM(:,:,3),coeffM(:,:,4),coeffM(:,:,5),coeffM(:,:,6),ref_lattice) +! ! Calculate the target quantities. All combinations of scalar products of +! ! the lattice vectors, b in XAX=b +! bvec(1)=dot_product(lattice(:,1),lattice(:,1)) +! bvec(2)=dot_product(lattice(:,2),lattice(:,2)) +! bvec(3)=dot_product(lattice(:,3),lattice(:,3)) +! bvec(4)=dot_product(lattice(:,1),lattice(:,2)) +! bvec(5)=dot_product(lattice(:,2),lattice(:,3)) +! bvec(6)=dot_product(lattice(:,3),lattice(:,1)) +! ! Do we have a guess for the strain? +! if ( present(guess) ) then +! eta0(1)=guess(1,1) +! eta0(2)=guess(2,2) +! eta0(3)=guess(3,3) +! eta0(4)=guess(1,2) +! eta0(5)=guess(1,3) +! eta0(6)=guess(2,3) +! else +! ! Start from something not completely crazy +! m0=matmul(ref_lattice,lo_invert3x3matrix(lattice)) +! eta0(1)=m0(1,1) +! eta0(2)=m0(2,2) +! eta0(3)=m0(3,3) +! eta0(4)=m0(1,2) +! eta0(5)=m0(1,3) +! eta0(6)=m0(2,3) +! endif +! eta1=0.0_flyt +! grad=0.0_flyt +! +! outloop: do outiter=1,nouter +! ! Get the error and gradient +! e0=0.0_flyt +! grad=0.0_flyt +! do i=1,6 +! f0=dot_product(eta0,matmul(coeffM(:,:,i),eta0))-bvec(i) +! e0=e0+f0**2 +! grad=grad+matmul( coeffM(:,:,i)+transpose(coeffM(:,:,i)),eta0)*2*f0 +! enddo +! +! ! Line search in this direction +! ctr=0 +! step=defstep +! inloop: do initer=1,ninner +! ! Guess new point +! eta1=eta0-step*grad +! ! Get error at new point +! e1=0.0_flyt +! do i=1,6 +! e1=e1+( dot_product(eta1,matmul(coeffM(:,:,i),eta1))-bvec(i) )**2 +! enddo +! +! ! See how it went: +! if ( e1 .lt. errtol ) then +! ! We are done +! eta0=eta1 +! e0=e1 +! exit outloop +! elseif ( e1 .lt. e0 ) then +! ! Good step, keep going +! eta0=eta1 +! e0=e1 +! step=step*1.3_flyt +! else +! ! Bad step, decrease and try again +! step=step*0.1_flyt +! ctr=ctr+1 +! endif +! +! ! If we have decreased too many times, start over with new gradient +! if ( ctr .ge. maxctr ) exit inloop +! enddo inloop +! enddo outloop +! +! ! Convert to a strain matrix +! strain(1,1)=eta0(1) +! strain(2,2)=eta0(2) +! strain(3,3)=eta0(3) +! strain(1,2)=eta0(4) +! strain(2,1)=eta0(4) +! strain(1,3)=eta0(5) +! strain(3,1)=eta0(5) +! strain(2,3)=eta0(6) +! strain(3,2)=eta0(6) +! ! Get the rotation +! m0=matmul(strain,ref_lattice) +! rotation=matmul(m0,lo_invert3x3matrix(lattice)) +! rotation=rotation/lo_determ(rotation)**(1.0_flyt/3.0_flyt) +! +! contains +! +! ! These are magic matrices that define the quadratic form. Thanks Mathematica. +! subroutine cfvars(cm1,cm2,cm3,cm4,cm5,cm6,basis) +! real(flyt), dimension(6,6), intent(out) :: cm1,cm2,cm3,cm4,cm5,cm6 +! real(flyt), dimension(3,3), intent(in) :: basis +! real(flyt) :: ax,ay,az,bx,by,bz,cx,cy,cz +! ax=basis(1,1) +! ay=basis(2,1) +! az=basis(3,1) +! bx=basis(1,2) +! by=basis(2,2) +! bz=basis(3,2) +! cx=basis(1,3) +! cy=basis(2,3) +! cz=basis(3,3) +! +! cm1(1,1)=ax**2 +! cm1(1,2)=0 +! cm1(1,3)=0 +! cm1(1,4)=2*ax*ay +! cm1(1,5)=2*ax*az +! cm1(1,6)=0 +! cm1(2,1)=0 +! cm1(2,2)=ay**2 +! cm1(2,3)=0 +! cm1(2,4)=2*ax*ay +! cm1(2,5)=0 +! cm1(2,6)=2*ay*az +! cm1(3,1)=0 +! cm1(3,2)=0 +! cm1(3,3)=az**2 +! cm1(3,4)=0 +! cm1(3,5)=2*ax*az +! cm1(3,6)=2*ay*az +! cm1(4,1)=0 +! cm1(4,2)=0 +! cm1(4,3)=0 +! cm1(4,4)=ax**2 + ay**2 +! cm1(4,5)=2*ay*az +! cm1(4,6)=2*ax*az +! cm1(5,1)=0 +! cm1(5,2)=0 +! cm1(5,3)=0 +! cm1(5,4)=0 +! cm1(5,5)=ax**2 + az**2 +! cm1(5,6)=2*ax*ay +! cm1(6,1)=0 +! cm1(6,2)=0 +! cm1(6,3)=0 +! cm1(6,4)=0 +! cm1(6,5)=0 +! cm1(6,6)=ay**2 + az**2 +! +! cm2(1,1)=bx**2 +! cm2(1,2)=0 +! cm2(1,3)=0 +! cm2(1,4)=2*bx*by +! cm2(1,5)=2*bx*bz +! cm2(1,6)=0 +! cm2(2,1)=0 +! cm2(2,2)=by**2 +! cm2(2,3)=0 +! cm2(2,4)=2*bx*by +! cm2(2,5)=0 +! cm2(2,6)=2*by*bz +! cm2(3,1)=0 +! cm2(3,2)=0 +! cm2(3,3)=bz**2 +! cm2(3,4)=0 +! cm2(3,5)=2*bx*bz +! cm2(3,6)=2*by*bz +! cm2(4,1)=0 +! cm2(4,2)=0 +! cm2(4,3)=0 +! cm2(4,4)=bx**2 + by**2 +! cm2(4,5)=2*by*bz +! cm2(4,6)=2*bx*bz +! cm2(5,1)=0 +! cm2(5,2)=0 +! cm2(5,3)=0 +! cm2(5,4)=0 +! cm2(5,5)=bx**2 + bz**2 +! cm2(5,6)=2*bx*by +! cm2(6,1)=0 +! cm2(6,2)=0 +! cm2(6,3)=0 +! cm2(6,4)=0 +! cm2(6,5)=0 +! cm2(6,6)=by**2 + bz**2 +! +! cm3(1,1)=cx**2 +! cm3(1,2)=0 +! cm3(1,3)=0 +! cm3(1,4)=2*cx*cy +! cm3(1,5)=2*cx*cz +! cm3(1,6)=0 +! cm3(2,1)=0 +! cm3(2,2)=cy**2 +! cm3(2,3)=0 +! cm3(2,4)=2*cx*cy +! cm3(2,5)=0 +! cm3(2,6)=2*cy*cz +! cm3(3,1)=0 +! cm3(3,2)=0 +! cm3(3,3)=cz**2 +! cm3(3,4)=0 +! cm3(3,5)=2*cx*cz +! cm3(3,6)=2*cy*cz +! cm3(4,1)=0 +! cm3(4,2)=0 +! cm3(4,3)=0 +! cm3(4,4)=cx**2 + cy**2 +! cm3(4,5)=2*cy*cz +! cm3(4,6)=2*cx*cz +! cm3(5,1)=0 +! cm3(5,2)=0 +! cm3(5,3)=0 +! cm3(5,4)=0 +! cm3(5,5)=cx**2 + cz**2 +! cm3(5,6)=2*cx*cy +! cm3(6,1)=0 +! cm3(6,2)=0 +! cm3(6,3)=0 +! cm3(6,4)=0 +! cm3(6,5)=0 +! cm3(6,6)=cy**2 + cz**2 +! +! cm4(1,1)=ax*bx +! cm4(1,2)=0 +! cm4(1,3)=0 +! cm4(1,4)=ay*bx + ax*by +! cm4(1,5)=az*bx + ax*bz +! cm4(1,6)=0 +! cm4(2,1)=0 +! cm4(2,2)=ay*by +! cm4(2,3)=0 +! cm4(2,4)=ay*bx + ax*by +! cm4(2,5)=0 +! cm4(2,6)=az*by + ay*bz +! cm4(3,1)=0 +! cm4(3,2)=0 +! cm4(3,3)=az*bz +! cm4(3,4)=0 +! cm4(3,5)=az*bx + ax*bz +! cm4(3,6)=az*by + ay*bz +! cm4(4,1)=0 +! cm4(4,2)=0 +! cm4(4,3)=0 +! cm4(4,4)=ax*bx + ay*by +! cm4(4,5)=az*by + ay*bz +! cm4(4,6)=az*bx + ax*bz +! cm4(5,1)=0 +! cm4(5,2)=0 +! cm4(5,3)=0 +! cm4(5,4)=0 +! cm4(5,5)=ax*bx + az*bz +! cm4(5,6)=ay*bx + ax*by +! cm4(6,1)=0 +! cm4(6,2)=0 +! cm4(6,3)=0 +! cm4(6,4)=0 +! cm4(6,5)=0 +! cm4(6,6)=ay*by + az*bz +! +! cm5(1,1)=bx*cx +! cm5(1,2)=0 +! cm5(1,3)=0 +! cm5(1,4)=by*cx + bx*cy +! cm5(1,5)=bz*cx + bx*cz +! cm5(1,6)=0 +! cm5(2,1)=0 +! cm5(2,2)=by*cy +! cm5(2,3)=0 +! cm5(2,4)=by*cx + bx*cy +! cm5(2,5)=0 +! cm5(2,6)=bz*cy + by*cz +! cm5(3,1)=0 +! cm5(3,2)=0 +! cm5(3,3)=bz*cz +! cm5(3,4)=0 +! cm5(3,5)=bz*cx + bx*cz +! cm5(3,6)=bz*cy + by*cz +! cm5(4,1)=0 +! cm5(4,2)=0 +! cm5(4,3)=0 +! cm5(4,4)=bx*cx + by*cy +! cm5(4,5)=bz*cy + by*cz +! cm5(4,6)=bz*cx + bx*cz +! cm5(5,1)=0 +! cm5(5,2)=0 +! cm5(5,3)=0 +! cm5(5,4)=0 +! cm5(5,5)=bx*cx + bz*cz +! cm5(5,6)=by*cx + bx*cy +! cm5(6,1)=0 +! cm5(6,2)=0 +! cm5(6,3)=0 +! cm5(6,4)=0 +! cm5(6,5)=0 +! cm5(6,6)=by*cy + bz*cz +! +! cm6(1,1)=ax*cx +! cm6(1,2)=0 +! cm6(1,3)=0 +! cm6(1,4)=ay*cx + ax*cy +! cm6(1,5)=az*cx + ax*cz +! cm6(1,6)=0 +! cm6(2,1)=0 +! cm6(2,2)=ay*cy +! cm6(2,3)=0 +! cm6(2,4)=ay*cx + ax*cy +! cm6(2,5)=0 +! cm6(2,6)=az*cy + ay*cz +! cm6(3,1)=0 +! cm6(3,2)=0 +! cm6(3,3)=az*cz +! cm6(3,4)=0 +! cm6(3,5)=az*cx + ax*cz +! cm6(3,6)=az*cy + ay*cz +! cm6(4,1)=0 +! cm6(4,2)=0 +! cm6(4,3)=0 +! cm6(4,4)=ax*cx + ay*cy +! cm6(4,5)=az*cy + ay*cz +! cm6(4,6)=az*cx + ax*cz +! cm6(5,1)=0 +! cm6(5,2)=0 +! cm6(5,3)=0 +! cm6(5,4)=0 +! cm6(5,5)=ax*cx + az*cz +! cm6(5,6)=ay*cx + ax*cy +! cm6(6,1)=0 +! cm6(6,2)=0 +! cm6(6,3)=0 +! cm6(6,4)=0 +! cm6(6,5)=0 +! cm6(6,6)=ay*cy + az*cz +!end subroutine +! +!end subroutine From d5eae40c20e956da793dc45d0ae673de7976dc2a Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Thu, 12 Dec 2024 14:52:46 +0100 Subject: [PATCH 02/11] average_structure | sanity check and graceful exit when no internal degrees of freedom are found --- src/average_structure/main.f90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/average_structure/main.f90 b/src/average_structure/main.f90 index 99d3acf4..0e04ab9e 100644 --- a/src/average_structure/main.f90 +++ b/src/average_structure/main.f90 @@ -44,6 +44,17 @@ program average_structure ! call sim%remove_force_and_center_of_mass_drift() call ic%generate(uc, verbosity=2, reference_positions=ur%r) + + ! report + print * + print *, '... No. of lattice degrees of freedom: ', tochar(ic%nx_lattice) + print *, '... No. of internal degrees of freedom: ', tochar(ic%nx_internal) + + ! Check if there is anything to average + if (ic%nx_internal .eq. 0) then + call lo_stop_gracefully(['*** No. of internal degrees of freedom is 0, no averaging needed'], 3) + end if + call avg_structure_nvt(uc, ss, ic, sim, nuc, nss, ur) !stop From b6e84fd6a84b8fc8ed199faec75db5c6ced77987 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 20:50:56 +0200 Subject: [PATCH 03/11] average_structure | use global optimization flags, fix Makefile deps, fix help text Co-Authored-By: Claude Opus 5.5 --- src/average_structure/Makefile | 5 ++--- src/average_structure/options.f90 | 6 +++--- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/average_structure/Makefile b/src/average_structure/Makefile index df6149ee..15a36083 100644 --- a/src/average_structure/Makefile +++ b/src/average_structure/Makefile @@ -9,7 +9,6 @@ LPATH = -L../../lib $(blaslapackLPATH) $(incLPATHhdf) $(incLPATHmpi) IPATH = -I../../inc/libolle -I../../inc/libflap $(blaslapackIPATH) $(incIPATHhdf) $(incIPATHmpi) LIBS = -lolle -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) -OPT = -O0 -fcheck=all -fbacktrace -fno-realloc-lhs -finit-real=nan -finit-derived -fmax-errors=10 F90 = $(FC) $(LPATH) $(IPATH) $(MODULE_FLAG) $(OBJECT_PATH) F90FLAGS = $(OPT) $(MODS) $(LIBS) @@ -19,14 +18,14 @@ $(PROG): $(OBJS) $(F90) $(LDFLAGS) -o $@ $(OBJS) $(LIBS) clean: - rm -f $(PROG) $(OBJS) *.mod + rm -f $(PROG) $(OBJS) $(OBJECT_PATH)*.mod .f90.o: $(F90) $(F90FLAGS) -c $< $(LIBS) $(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o $(OBJECT_PATH)type_irredcell.o $(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@ -$(OBJECT_PATH)type_irredcell.o: +$(OBJECT_PATH)type_irredcell.o: type_irredcell.f90 type_irredcell_aux.f90 $(F90) $(F90FLAGS) -c type_irredcell.f90 $(LIBS) -o $@ $(OBJECT_PATH)options.o: $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ diff --git a/src/average_structure/options.f90 b/src/average_structure/options.f90 index 294feb1a..bb0e8d8a 100644 --- a/src/average_structure/options.f90 +++ b/src/average_structure/options.f90 @@ -35,14 +35,14 @@ subroutine parse(opts) version=lo_version, & license=lo_licence, & help='Usage: ', & - description='Utility to fit equations of state to volume-energy data.', & - examples=["eosfit"], & + description='Compute the symmetry-respecting average structure from an MD simulation.', & + examples=["average_structure"], & epilog=new_line('a')//"...") cli_manpage cli_verbose call cli%add(switch='--stride', switch_ab='-s', & - help='Number of dimenions for the fit', & + help='Use every N configuration instead of all.', & required=.false., act='store', def='1', error=lo_status) if (lo_status .ne. 0) stop From 16a8163b347a8a643983ee733bc50525e63dfdbd Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:21:54 +0200 Subject: [PATCH 04/11] average_structure | add test for symmetry-projected average structure Co-Authored-By: Claude Opus 5.5 --- tests/Makefile | 1 + tests/average_structure/Makefile | 6 ++ .../test_average_structure.py | 59 +++++++++++++++++++ .../tetragonal_ab_average/infile.forces | 4 ++ .../tetragonal_ab_average/infile.meta | 4 ++ .../tetragonal_ab_average/infile.positions | 4 ++ .../tetragonal_ab_average/infile.refposcar | 10 ++++ .../tetragonal_ab_average/infile.ssposcar | 10 ++++ .../tetragonal_ab_average/infile.stat | 2 + .../tetragonal_ab_average/infile.ucposcar | 10 ++++ tests/make_all_testfiles.sh | 1 + 11 files changed, 111 insertions(+) create mode 100644 tests/average_structure/Makefile create mode 100644 tests/average_structure/test_average_structure.py create mode 100644 tests/average_structure/tetragonal_ab_average/infile.forces create mode 100644 tests/average_structure/tetragonal_ab_average/infile.meta create mode 100644 tests/average_structure/tetragonal_ab_average/infile.positions create mode 100644 tests/average_structure/tetragonal_ab_average/infile.refposcar create mode 100644 tests/average_structure/tetragonal_ab_average/infile.ssposcar create mode 100644 tests/average_structure/tetragonal_ab_average/infile.stat create mode 100644 tests/average_structure/tetragonal_ab_average/infile.ucposcar diff --git a/tests/Makefile b/tests/Makefile index 22a95542..210743b8 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -9,6 +9,7 @@ all: clean testfiles clean: rm -f */outfile* + rm -f */*/outfile* rm -rf */__pycache__/ rm -rf *.pytest_cache/ rm -f */contcar* diff --git a/tests/average_structure/Makefile b/tests/average_structure/Makefile new file mode 100644 index 00000000..5016e0fd --- /dev/null +++ b/tests/average_structure/Makefile @@ -0,0 +1,6 @@ +cases = tetragonal_ab_average + +testfiles: + for case in $(cases); do \ + (cd $$case && average_structure > outfile.log 2>&1; echo $$? > outfile.exitcode); \ + done diff --git a/tests/average_structure/test_average_structure.py b/tests/average_structure/test_average_structure.py new file mode 100644 index 00000000..c6209244 --- /dev/null +++ b/tests/average_structure/test_average_structure.py @@ -0,0 +1,59 @@ +import re + +import numpy as np +from pathlib import Path + +parent = Path(__file__).parent + + +def _exitcode(case): + return int((parent / case / "outfile.exitcode").read_text()) + + +def _n_internal_dof(case): + log = (parent / case / "outfile.log").read_text() + match = re.search(r"No\. of internal degrees of freedom:\s*(\d+)", log) + assert match, f"no DOF count in {parent / case / 'outfile.log'}" + return int(match.group(1)) + + +def _read_poscar(file): + lines = Path(file).read_text().splitlines() + lattice = np.array([line.split() for line in lines[2:5]], dtype=float) + lattice *= float(lines[1]) + species = lines[5].split() + counts = [int(n) for n in lines[6].split()] + symbols = [s for s, n in zip(species, counts) for _ in range(n)] + positions = np.array([line.split()[:3] for line in lines[8 : 8 + sum(counts)]], dtype=float) + return lattice, symbols, positions + + +def test_average_is_symmetry_projected(): + """ + Scenario: the average structure is the symmetry-projected average of the trajectory + Given a tetragonal AB crystal (P4mm, A = Ga, B = N) with A at (0,0,0) and B at (1/2,1/2,0.600) + And an MD trajectory of two frames: + frame 1 = B at z = 0.580 + frame 2 = B at z = 0.640, and A additionally shifted by (+0.010, 0, 0) + When the average structure is computed + Then 1 internal degree of freedom is reported + And outfile.ucposcar has A at (0,0,-0.005) and B at (1/2,1/2,0.605) + And the lattice is unchanged + """ + case = "tetragonal_ab_average" + assert _exitcode(case) == 0 + assert _n_internal_dof(case) == 1 + + lattice_ref, _, _ = _read_poscar(parent / case / "infile.ucposcar") + lattice, _, positions = _read_poscar(parent / case / "outfile.ucposcar") + np.testing.assert_allclose(lattice, lattice_ref, atol=1e-8) + + positions_expected = np.array([[0.0, 0.0, -0.005], [0.5, 0.5, 0.605]]) + # compare modulo lattice translations + diff = positions - positions_expected + diff -= np.round(diff) + np.testing.assert_allclose(diff, 0.0, atol=1e-8) + + +if __name__ == "__main__": + test_average_is_symmetry_projected() diff --git a/tests/average_structure/tetragonal_ab_average/infile.forces b/tests/average_structure/tetragonal_ab_average/infile.forces new file mode 100644 index 00000000..b817a0ab --- /dev/null +++ b/tests/average_structure/tetragonal_ab_average/infile.forces @@ -0,0 +1,4 @@ + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 diff --git a/tests/average_structure/tetragonal_ab_average/infile.meta b/tests/average_structure/tetragonal_ab_average/infile.meta new file mode 100644 index 00000000..1e50194b --- /dev/null +++ b/tests/average_structure/tetragonal_ab_average/infile.meta @@ -0,0 +1,4 @@ + 2 # N atoms + 2 # N timesteps + 1.0 # timestep in fs + 300.0 # temperature in K diff --git a/tests/average_structure/tetragonal_ab_average/infile.positions b/tests/average_structure/tetragonal_ab_average/infile.positions new file mode 100644 index 00000000..de652968 --- /dev/null +++ b/tests/average_structure/tetragonal_ab_average/infile.positions @@ -0,0 +1,4 @@ + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 5.000000000000000e-01 5.000000000000000e-01 5.800000000000000e-01 + 1.000000000000000e-02 0.000000000000000e+00 0.000000000000000e+00 + 5.000000000000000e-01 5.000000000000000e-01 6.400000000000000e-01 diff --git a/tests/average_structure/tetragonal_ab_average/infile.refposcar b/tests/average_structure/tetragonal_ab_average/infile.refposcar new file mode 100644 index 00000000..69d4144f --- /dev/null +++ b/tests/average_structure/tetragonal_ab_average/infile.refposcar @@ -0,0 +1,10 @@ + Ga N + 1.0000000000000000 + 4.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 4.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 4.0000000000000000 + Ga N + 1 1 +Direct + 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.5000000000000000 0.5000000000000000 0.6000000000000000 diff --git a/tests/average_structure/tetragonal_ab_average/infile.ssposcar b/tests/average_structure/tetragonal_ab_average/infile.ssposcar new file mode 100644 index 00000000..69d4144f --- /dev/null +++ b/tests/average_structure/tetragonal_ab_average/infile.ssposcar @@ -0,0 +1,10 @@ + Ga N + 1.0000000000000000 + 4.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 4.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 4.0000000000000000 + Ga N + 1 1 +Direct + 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.5000000000000000 0.5000000000000000 0.6000000000000000 diff --git a/tests/average_structure/tetragonal_ab_average/infile.stat b/tests/average_structure/tetragonal_ab_average/infile.stat new file mode 100644 index 00000000..9c500354 --- /dev/null +++ b/tests/average_structure/tetragonal_ab_average/infile.stat @@ -0,0 +1,2 @@ + 1 0.000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 2 1.000 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 diff --git a/tests/average_structure/tetragonal_ab_average/infile.ucposcar b/tests/average_structure/tetragonal_ab_average/infile.ucposcar new file mode 100644 index 00000000..69d4144f --- /dev/null +++ b/tests/average_structure/tetragonal_ab_average/infile.ucposcar @@ -0,0 +1,10 @@ + Ga N + 1.0000000000000000 + 4.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.0000000000000000 4.0000000000000000 0.0000000000000000 + 0.0000000000000000 0.0000000000000000 4.0000000000000000 + Ga N + 1 1 +Direct + 0.0000000000000000 0.0000000000000000 0.0000000000000000 + 0.5000000000000000 0.5000000000000000 0.6000000000000000 diff --git a/tests/make_all_testfiles.sh b/tests/make_all_testfiles.sh index 3292df39..2986c851 100644 --- a/tests/make_all_testfiles.sh +++ b/tests/make_all_testfiles.sh @@ -1,5 +1,6 @@ folders="anharmonic_free_energy/ atomic_distribution/ +average_structure/ canonical_configuration/ crystal_structure_info/ dump_dynamical_matrices/ From 6454e16aaf752d77986db60ef3b6a19b3801649f Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:27:13 +0200 Subject: [PATCH 05/11] average_structure | make infile.refposcar optional and check it against infile.ucposcar Co-Authored-By: Claude Opus 5.5 --- src/average_structure/main.f90 | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/average_structure/main.f90 b/src/average_structure/main.f90 index 0e04ab9e..3ede38ea 100644 --- a/src/average_structure/main.f90 +++ b/src/average_structure/main.f90 @@ -34,8 +34,20 @@ program average_structure call ss%readfromfile('infile.ssposcar') call ss%classify('supercell', uc) - ! Reference positions to determine displacement from - call ur%readfromfile('infile.refposcar') + ! Reference positions to determine displacement from, default to the unitcell + if (lo_does_file_exist('infile.refposcar')) then + call ur%readfromfile('infile.refposcar') + if (ur%na .ne. uc%na) then + call lo_stop_gracefully(['infile.refposcar and infile.ucposcar have different number of atoms'], & + lo_exitcode_param, __FILE__, __LINE__) + end if + if (any(ur%atomic_number .ne. uc%atomic_number)) then + call lo_stop_gracefully(['infile.refposcar and infile.ucposcar have different species or atom order'], & + lo_exitcode_param, __FILE__, __LINE__) + end if + else + ur = uc + end if ! Read the MD simulation call sim%read_from_file(verbosity=2, stride=opts%stride) From 6492e5a096dbd48d8972936a63413d3c9a5089c3 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:27:26 +0200 Subject: [PATCH 06/11] average_structure | check that --stride is a positive integer Co-Authored-By: Claude Opus 5.5 --- src/average_structure/options.f90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/average_structure/options.f90 b/src/average_structure/options.f90 index bb0e8d8a..7be1714b 100644 --- a/src/average_structure/options.f90 +++ b/src/average_structure/options.f90 @@ -1,6 +1,7 @@ #include "precompilerdefinitions" module options -use konstanter, only: flyt, lo_status, lo_author, lo_version, lo_licence, lo_hugeint +use konstanter, only: flyt, lo_status, lo_author, lo_version, lo_licence, lo_hugeint, lo_exitcode_param +use gottochblandat, only: lo_stop_gracefully use flap, only: command_line_interface implicit none @@ -64,6 +65,10 @@ subroutine parse(opts) ! Parse the rest call cli%get(switch='--stride', val=opts%stride, error=lo_status) + if (lo_status .ne. 0) stop + if (opts%stride .lt. 1) then + call lo_stop_gracefully(['--stride has to be a positive integer'], lo_exitcode_param, __FILE__, __LINE__) + end if end subroutine end module From ebfac8639558b2813b78cbf4275a97368d11d195 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:27:45 +0200 Subject: [PATCH 07/11] average_structure | write the unchanged structure when there are no internal degrees of freedom Co-Authored-By: Claude Opus 5.5 --- src/average_structure/main.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/average_structure/main.f90 b/src/average_structure/main.f90 index 3ede38ea..8d537bd9 100644 --- a/src/average_structure/main.f90 +++ b/src/average_structure/main.f90 @@ -62,13 +62,15 @@ program average_structure print *, '... No. of lattice degrees of freedom: ', tochar(ic%nx_lattice) print *, '... No. of internal degrees of freedom: ', tochar(ic%nx_internal) - ! Check if there is anything to average + ! Check if there is anything to average, otherwise the structure is fixed by symmetry if (ic%nx_internal .eq. 0) then - call lo_stop_gracefully(['*** No. of internal degrees of freedom is 0, no averaging needed'], 3) + write (*, *) '... No. of internal degrees of freedom is 0, no averaging needed' + nuc = uc + nss = ss + else + call avg_structure_nvt(uc, ss, ic, sim, nuc, nss, ur) end if - call avg_structure_nvt(uc, ss, ic, sim, nuc, nss, ur) - !stop ! call sim%remove_force_and_center_of_mass_drift() ! ! Remove arbitrary rotations From 08a98b6700d4042f0a1aea74650690a7d651f73d Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:27:58 +0200 Subject: [PATCH 08/11] average_structure | rebuild main.o and options.o when their sources change Co-Authored-By: Claude Opus 5.5 --- src/average_structure/Makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/average_structure/Makefile b/src/average_structure/Makefile index 15a36083..69eeef48 100644 --- a/src/average_structure/Makefile +++ b/src/average_structure/Makefile @@ -23,10 +23,10 @@ clean: .f90.o: $(F90) $(F90FLAGS) -c $< $(LIBS) -$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o $(OBJECT_PATH)type_irredcell.o +$(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o $(OBJECT_PATH)type_irredcell.o $(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@ $(OBJECT_PATH)type_irredcell.o: type_irredcell.f90 type_irredcell_aux.f90 $(F90) $(F90FLAGS) -c type_irredcell.f90 $(LIBS) -o $@ -$(OBJECT_PATH)options.o: +$(OBJECT_PATH)options.o: options.f90 $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ From 164370c97c86c664fae4791afd0456fb0ac16217 Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:28:18 +0200 Subject: [PATCH 09/11] average_structure | remove unused npt option, warn that infile.lattice is ignored Co-Authored-By: Claude Opus 5.5 --- src/average_structure/main.f90 | 5 ++++- src/average_structure/options.f90 | 4 +--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/average_structure/main.f90 b/src/average_structure/main.f90 index 8d537bd9..a47c18cf 100644 --- a/src/average_structure/main.f90 +++ b/src/average_structure/main.f90 @@ -49,8 +49,11 @@ program average_structure ur = uc end if - ! Read the MD simulation + ! Read the MD simulation, only fixed-cell simulations are supported call sim%read_from_file(verbosity=2, stride=opts%stride) + if (lo_does_file_exist('infile.lattice')) then + write (*, *) '... WARNING: infile.lattice is ignored, the lattice is taken from infile.ucposcar' + end if ! lo_deallocate(sim%lattice) ! Remove center of mass drift ! call sim%remove_force_and_center_of_mass_drift() diff --git a/src/average_structure/options.f90 b/src/average_structure/options.f90 index 7be1714b..3330319f 100644 --- a/src/average_structure/options.f90 +++ b/src/average_structure/options.f90 @@ -13,8 +13,6 @@ module options integer :: verbosity = -lo_hugeint !> use every N timesteps integer :: stride = -lo_hugeint - !> is it a constant pressure simulation? - logical :: npt = .false. contains procedure :: parse end type @@ -36,7 +34,7 @@ subroutine parse(opts) version=lo_version, & license=lo_licence, & help='Usage: ', & - description='Compute the symmetry-respecting average structure from an MD simulation.', & + description='Compute the symmetry-respecting average structure from an MD simulation with fixed cell (NVT).', & examples=["average_structure"], & epilog=new_line('a')//"...") From f684624a8ffe56bb078b19864c6ccac150c684cb Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:28:39 +0200 Subject: [PATCH 10/11] average_structure | add documentation Co-Authored-By: Claude Opus 5.5 --- docs/program/README.md | 2 ++ docs/program/average_structure.md | 48 +++++++++++++++++++++++++++++++ mkdocs.yml | 1 + src/average_structure/manual.md | 23 +++++++++++++++ 4 files changed, 74 insertions(+) create mode 100644 docs/program/average_structure.md create mode 100644 src/average_structure/manual.md diff --git a/docs/program/README.md b/docs/program/README.md index fff2eb91..060ad66a 100644 --- a/docs/program/README.md +++ b/docs/program/README.md @@ -20,6 +20,8 @@ [`samples_from_md`: Pick samples from an MD simulation in a clever way.](samples_from_md.md) +[`average_structure`: Compute the symmetry-respecting average structure from an MD simulation.](average_structure.md) + [`dump_dynamical_matrices`: Write dynamical matrices on q-point grid.](dump_dynamical_matrices.md) [`crystal_structure_info`: Report which crystal structure and spacegroup TDEP sees.](crystal_structure_info.md) diff --git a/docs/program/average_structure.md b/docs/program/average_structure.md new file mode 100644 index 00000000..d8f2e88b --- /dev/null +++ b/docs/program/average_structure.md @@ -0,0 +1,48 @@ + +### Short description + +Compute the symmetry-respecting average structure from an MD simulation with fixed cell (NVT). + +### Command line options: + + + + +Optional switches: + +* `--stride value`, `-s value` + default value 1 + Use every N configuration instead of all. + +* `--help`, `-h` + Print this help message + +* `--version`, `-v` + Print version +### Examples + +`average_structure` + +### Longer summary + +This code computes the average structure from a molecular dynamics simulation such that the result has the symmetry of the reference structure. The internal degrees of freedom, i.e. the atomic displacements that are allowed by the space group of `infile.ucposcar`, are determined first. The net translation of the cell is removed by requiring that the displacements sum to zero. + +For each timestep, the displacements of all atoms in the supercell from the reference positions are projected onto these internal degrees of freedom. The projections are averaged over the trajectory and the averaged structure is constructed from them. Displacements that break the symmetry are projected out, and the averaged structure keeps the space group of the reference. The lattice is not changed, so only simulations with fixed cell (NVT) are supported. + +If the structure has no internal degrees of freedom, i.e. all positions are fixed by symmetry, the input structure is written unchanged. + +### Input files + +* [infile.ucposcar](../files.md#infile.ucposcar), determines the symmetry +* [infile.ssposcar](../files.md#infile.ssposcar) +* [infile.meta](../files.md#infile.meta) +* [infile.stat](../files.md#infile.stat) +* [infile.positions](../files.md#infile.positions) +* [infile.forces](../files.md#infile.forces) +* `infile.refposcar` (optional), the reference positions the displacements are measured from. Same format, number of atoms, species and atom order as `infile.ucposcar`. Defaults to the positions in `infile.ucposcar`. + +### Output files + +* `outfile.ucposcar`, the averaged unit cell +* `outfile.ssposcar`, the averaged supercell, consistent with `infile.ssposcar` +* `outfile.internal`, the time in fs (first column) and the value of each internal degree of freedom (remaining columns) for every timestep diff --git a/mkdocs.yml b/mkdocs.yml index c477875c..9f49e316 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -40,6 +40,7 @@ nav: - Atomic distribution: program/atomic_distribution.md - Pack simulation: program/pack_simulation.md - Samples from MD: program/samples_from_md.md + - Average structure: program/average_structure.md - Dump dynamical matrices: program/dump_dynamical_matrices.md - Crystal structure info: program/crystal_structure_info.md - Refine structure: program/refine_structure.md diff --git a/src/average_structure/manual.md b/src/average_structure/manual.md new file mode 100644 index 00000000..e6495747 --- /dev/null +++ b/src/average_structure/manual.md @@ -0,0 +1,23 @@ +### Longer summary + +This code computes the average structure from a molecular dynamics simulation such that the result has the symmetry of the reference structure. The internal degrees of freedom, i.e. the atomic displacements that are allowed by the space group of `infile.ucposcar`, are determined first. The net translation of the cell is removed by requiring that the displacements sum to zero. + +For each timestep, the displacements of all atoms in the supercell from the reference positions are projected onto these internal degrees of freedom. The projections are averaged over the trajectory and the averaged structure is constructed from them. Displacements that break the symmetry are projected out, and the averaged structure keeps the space group of the reference. The lattice is not changed, so only simulations with fixed cell (NVT) are supported. + +If the structure has no internal degrees of freedom, i.e. all positions are fixed by symmetry, the input structure is written unchanged. + +### Input files + +* [infile.ucposcar](../files.md#infile.ucposcar), determines the symmetry +* [infile.ssposcar](../files.md#infile.ssposcar) +* [infile.meta](../files.md#infile.meta) +* [infile.stat](../files.md#infile.stat) +* [infile.positions](../files.md#infile.positions) +* [infile.forces](../files.md#infile.forces) +* `infile.refposcar` (optional), the reference positions the displacements are measured from. Same format, number of atoms, species and atom order as `infile.ucposcar`. Defaults to the positions in `infile.ucposcar`. + +### Output files + +* `outfile.ucposcar`, the averaged unit cell +* `outfile.ssposcar`, the averaged supercell, consistent with `infile.ssposcar` +* `outfile.internal`, the time in fs (first column) and the value of each internal degree of freedom (remaining columns) for every timestep From 046a9d776d152745dfa51aded8a459776b604a6b Mon Sep 17 00:00:00 2001 From: Florian Knoop Date: Wed, 7 Oct 2026 21:34:26 +0200 Subject: [PATCH 11/11] average_structure | link libolle statically like the other codes Co-Authored-By: Claude Opus 5.5 --- src/average_structure/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/average_structure/Makefile b/src/average_structure/Makefile index 69eeef48..1438f944 100644 --- a/src/average_structure/Makefile +++ b/src/average_structure/Makefile @@ -7,7 +7,7 @@ OBJS = $(OBJECT_PATH)main.o $(OBJECT_PATH)options.o $(OBJECT_PATH)type_irredcell LPATH = -L../../lib $(blaslapackLPATH) $(incLPATHhdf) $(incLPATHmpi) IPATH = -I../../inc/libolle -I../../inc/libflap $(blaslapackIPATH) $(incIPATHhdf) $(incIPATHmpi) -LIBS = -lolle -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) +LIBS = ../../lib/libolle.a -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) F90 = $(FC) $(LPATH) $(IPATH) $(MODULE_FLAG) $(OBJECT_PATH) F90FLAGS = $(OPT) $(MODS) $(LIBS)