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/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/Makefile b/src/average_structure/Makefile new file mode 100644 index 00000000..1438f944 --- /dev/null +++ b/src/average_structure/Makefile @@ -0,0 +1,32 @@ +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 = ../../lib/libolle.a -lflap $(blaslapackLIBS) $(incLIBShdf) $(incLIBSmpi) + +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) $(OBJECT_PATH)*.mod + +.f90.o: + $(F90) $(F90FLAGS) -c $< $(LIBS) + +$(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: options.f90 + $(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..a47c18cf --- /dev/null +++ b/src/average_structure/main.f90 @@ -0,0 +1,273 @@ +#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, 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, 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() + + 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, otherwise the structure is fixed by symmetry + if (ic%nx_internal .eq. 0) then + 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 + +!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/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 diff --git a/src/average_structure/options.f90 b/src/average_structure/options.f90 new file mode 100644 index 00000000..3330319f --- /dev/null +++ b/src/average_structure/options.f90 @@ -0,0 +1,72 @@ +#include "precompilerdefinitions" +module options +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 + +private +public :: lo_opts + +type lo_opts + !> how much to talk + integer :: verbosity = -lo_hugeint + !> use every N timesteps + integer :: stride = -lo_hugeint +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='Compute the symmetry-respecting average structure from an MD simulation with fixed cell (NVT).', & + examples=["average_structure"], & + epilog=new_line('a')//"...") + + cli_manpage + cli_verbose + call cli%add(switch='--stride', switch_ab='-s', & + help='Use every N configuration instead of all.', & + 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) + 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 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 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/