diff --git a/src/canonical_configuration/main.f90 b/src/canonical_configuration/main.f90 index c5aaca52..3b4845e1 100644 --- a/src/canonical_configuration/main.f90 +++ b/src/canonical_configuration/main.f90 @@ -210,6 +210,12 @@ program canonical_configuration case (5) ! Siesta output fname = 'siesta_conf'//tochar(iconf, 4) call p%writetofile(trim(fname), opts%output_format, write_velocities=.true.) + case (6) ! QE output + fname = 'qe_conf'//tochar(iconf, 4) + call p%writetofile(trim(fname), opts%output_format, write_velocities=.false.) + case (7) ! QE output + fname = 'parsec_conf'//tochar(iconf, 4) + call p%writetofile(trim(fname), opts%output_format, write_velocities=.false.) end select ! just measure some stuff, for no good reason. diff --git a/src/canonical_configuration/options.f90 b/src/canonical_configuration/options.f90 index db68e36c..cbe3d257 100644 --- a/src/canonical_configuration/options.f90 +++ b/src/canonical_configuration/options.f90 @@ -66,9 +66,9 @@ subroutine parse(opts) help_markdown="That is, use \( \sqrt{\frac{\hbar (2n+1) }{2 m \omega}} \) as the mean normal mode amplitudes instead of the classical \( \frac{1}{\omega}\sqrt{\frac{k_BT}{m}} \)", & required=.false., act='store_true', def='.false.', error=lo_status) if (lo_status .ne. 0) stop - call cli%add(switch='--output_format', switch_ab='-of', & - help='Selects output format. 1 is VASP, 2 is Abinit, 4 is FHI-Aims, 5 is Siesta. Default 1.', & - required=.false., act='store', def='1', choices='1,2,4,5', error=lo_status) + call cli%add(switch='--output_format', switch_ab='-of', hidden=.false., & + help='Output format. 1 is VASP, 2 Abinit, 4 FHI-Aims, 5 Siesta, 6 QE and 7 Parsec', & + required=.false., act='store', def='1', choices='1,2,4,5,6,7', error=lo_status) if (lo_status .ne. 0) stop call cli%add(switch='--mindist', & help='What is the smallest distance between two atoms allowed, in units of the nearest neighbour distance.', & diff --git a/src/dump_dynamical_matrices/Makefile b/src/dump_dynamical_matrices/Makefile index 7c6d6b2a..47065b30 100644 --- a/src/dump_dynamical_matrices/Makefile +++ b/src/dump_dynamical_matrices/Makefile @@ -23,7 +23,7 @@ clean: rm -f $(PROG) $(OBJS) *.mod $(OBJECT_PATH)*.mod -$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o +$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o ../../lib/libolle.a $(F90) $(F90FLAGS) -c main.f90 -o $@ $(OBJECT_PATH)options.o: $(F90) $(F90FLAGS) -c options.f90 -o $@ diff --git a/src/dump_dynamical_matrices/main.f90 b/src/dump_dynamical_matrices/main.f90 index 841ec7f5..1402dfae 100644 --- a/src/dump_dynamical_matrices/main.f90 +++ b/src/dump_dynamical_matrices/main.f90 @@ -78,6 +78,8 @@ program dump_dynamical_matrices call fc%write_to_anaddb(uc, opts%qgrid, mw, mem) +call fc%write_dynmat_to_qe(uc,opts%qgrid,mw,mem) + !else ! ! automagically generate a grid. Probably a bad idea, since my version of gridgeneration and ! ! Abinits might be different. diff --git a/src/extract_forceconstants/Makefile b/src/extract_forceconstants/Makefile index 55290768..1e14493f 100644 --- a/src/extract_forceconstants/Makefile +++ b/src/extract_forceconstants/Makefile @@ -26,7 +26,7 @@ clean: .f90.o: $(F90) $(F90FLAGS) -c $< $(LIBS) -$(OBJECT_PATH)main.o: $(OBJECT_PATH)options.o +$(OBJECT_PATH)main.o: main.f90 $(OBJECT_PATH)options.o $(F90) $(F90FLAGS) -c main.f90 $(LIBS) -o $@ -$(OBJECT_PATH)options.o: +$(OBJECT_PATH)options.o: options.f90 $(F90) $(F90FLAGS) -c options.f90 $(LIBS) -o $@ diff --git a/src/extract_forceconstants/main.f90 b/src/extract_forceconstants/main.f90 index 60443d8c..e94a2bb3 100644 --- a/src/extract_forceconstants/main.f90 +++ b/src/extract_forceconstants/main.f90 @@ -348,6 +348,11 @@ program extract_forceconstants do i = 1, 6 write (*, "(6(3X,F15.5))") lo_chop((fc2%elastic_constants_voigt(:, i) + fc2%elastic_constants_voigt_longrange(:, i))*lo_pressure_HartreeBohr_to_GPa, lo_tol) end do + + if ( opts%dumpepw ) then + ! If desired, write IFCs in EPW format + call fc2%write_to_qe(uc,mw,mem) + endif end if end if if (map%have_fc_triplet) then diff --git a/src/extract_forceconstants/options.f90 b/src/extract_forceconstants/options.f90 index 7d1f1ca7..043ccbce 100644 --- a/src/extract_forceconstants/options.f90 +++ b/src/extract_forceconstants/options.f90 @@ -68,7 +68,8 @@ module options logical :: fake_dielectric = .false. ! developer mode logical :: devmode = .false. - + ! Dump IFCs in EPW format + logical :: dumpepw = .false. contains procedure :: parse end type @@ -248,6 +249,9 @@ subroutine parse(opts) call cli%add(switch='--developermode', switch_ab='-dev', hidden=.true., help='dev. mode', & required=.false., act='store_true', def='.false.', error=lo_status) if (lo_status .ne. 0) stop + call cli%add(switch='--output_epw', switch_ab='-epw', hidden=.false., help='Write the (second order) interatomic force constants in XML format useable by EPW.', & + required=.false., act='store_true', def='.false.', error=lo_status) + if (lo_status .ne. 0) stop ! actually parse it call cli%parse(error=lo_status) @@ -306,6 +310,7 @@ subroutine parse(opts) call cli%get(switch='--temperature', val=opts%temperature) call cli%get(switch='--fakediel', val=opts%fake_dielectric) call cli%get(switch='--developermode', val=opts%devmode) + call cli%get(switch='--output_epw', val=opts%dumpepw) if (lo_status .ne. 0) stop diff --git a/src/generate_structure/main.f90 b/src/generate_structure/main.f90 index e3825c13..d004db11 100644 --- a/src/generate_structure/main.f90 +++ b/src/generate_structure/main.f90 @@ -113,6 +113,13 @@ program generate_structure ss%latticevectors = ss%latticevectors ss%inv_latticevectors = ss%inv_latticevectors call ss%writetofile('outfile.supercell_lammps_ipi', 3, transformationmatrix=tm) + case (6) ! QE + call ss%writetofile('outfile.supercell_qe', opts%outputformat) + write (*, *) '... wrote supercell in Quantum ESPRESSO format' + case (7) ! Parsec + call uc%writetofile('outfile.unitcell_parsec', opts%outputformat) + call ss%writetofile('outfile.supercell_parsec', opts%outputformat) + write (*, *) '... wrote supercell in PARSEC format' end select end if end block getsupercell diff --git a/src/generate_structure/options.f90 b/src/generate_structure/options.f90 index da5e6e3f..f6c2beb7 100644 --- a/src/generate_structure/options.f90 +++ b/src/generate_structure/options.f90 @@ -64,8 +64,8 @@ subroutine parse(opts) required=.false., act='store', def='5.0', error=lo_status) if (lo_status .ne. 0) stop call cli%add(switch='--output_format', switch_ab='-of', hidden=.false., & - help='Output format. 1 is VASP, 2 Abinit, 4 FHI-Aims and 5 Siesta', & - required=.false., act='store', def='1', choices='1,2,4,5', error=lo_status) + help='Output format. 1 is VASP, 2 Abinit, 4 FHI-Aims, 5 Siesta, 6 QE and 7 Parsec', & + required=.false., act='store', def='1', choices='1,2,4,5,6,7', error=lo_status) if (lo_status .ne. 0) stop call cli%add(switch='--thirdorder_cutoff', switch_ab='-rc3', hidden=.true., & help='Cutoff for the third order force constants', & diff --git a/src/libolle/type_crystalstructure_io.f90 b/src/libolle/type_crystalstructure_io.f90 index 771a61f9..465103ba 100644 --- a/src/libolle/type_crystalstructure_io.f90 +++ b/src/libolle/type_crystalstructure_io.f90 @@ -376,8 +376,12 @@ module subroutine writetofile(p,filename,output_format,write_velocities,transfor end if case(4) ! FHI Aims call writetofile_aims(p,filename,write_velocities) - case(5) ! Siesta + case(5) ! Siesta call writetofile_siesta(p,filename,write_velocities) + case(6) ! QE + call writetofile_qe(p,filename,write_velocities) + case(7) ! Parsec + call writetofile_parsec(p,filename,write_velocities) case default call lo_stop_gracefully(['Unknown output format: '//tochar(output_format)],lo_exitcode_io,__FILE__,__LINE__) end select @@ -726,7 +730,7 @@ subroutine writetofile_siesta(p,filename,write_velocities) ! I don't know how to write velocities in Siesta. Someone should tell me that. ! MJV: fixed at least partly, by putting them in XV files (pos, veloc). ! not sure you can add them to the main fdf input files as well... - ! NB: there is a flag to use the XV file MD.UseSaveXV .true. but you have to fill + ! NB: there is a flag to use the XV file MD.UseSaveXV .true. but you have to fill ! the coordinates block anyway, so might as well fill both here. ! Main advantage is that you can add the velocities in XV, but for configuration ! generation it's not important, you just get forces out. @@ -785,7 +789,7 @@ subroutine writetofile_siesta(p,filename,write_velocities) write(u,*) "XC.functional GGA" write(u,*) "XC.authors PBE" write(u,*) - + write(u,*) "%block kgrid_Monkhorst_Pack" write(u,*) " 1 0 0 0." write(u,*) " 0 1 0 0." @@ -803,25 +807,114 @@ subroutine writetofile_siesta(p,filename,write_velocities) if ( filename .ne. 'stdout' ) close(u) -! add a XV format file as well: Uses bohr for the cartesian positions. -! what are the units for the velocity? -! -! latice as 3x3 block, with additional 3x3 block of 0s after -! natom -! itype zatom xcoord(3) veloc(3) -! -! NB: I _think_ they use Angstr for the positions, needs to be checked. -! for the velocities, no idea what they expect... + ! add a XV format file as well: Uses bohr for the cartesian positions. + ! what are the units for the velocity? + ! + ! latice as 3x3 block, with additional 3x3 block of 0s after + ! natom + ! itype zatom xcoord(3) veloc(3) + ! + ! NB: I _think_ they use Angstr for the positions, needs to be checked. + ! for the velocities, no idea what they expect... u=open_file('out',filename//".XV") write(u,'(3E20.10,2x,3E20.10)') lo_choplarge(p%latticevectors(:,1),lo_sqtol), 0.000000000, 0.000000000, 0.000000000 write(u,'(3E20.10,2x,3E20.10)') lo_choplarge(p%latticevectors(:,2),lo_sqtol), 0.000000000, 0.000000000, 0.000000000 write(u,'(3E20.10,2x,3E20.10)') lo_choplarge(p%latticevectors(:,3),lo_sqtol), 0.000000000, 0.000000000, 0.000000000 - write(u,*) tochar(p%na) + write(u,*) tochar(p%na) opf="(1X,2(I5),2X,2(3F18.9,' '))" do i=1,p%na write(u,opf) p%species(i), symbol_to_z(p%atomic_symbol(p%species(i))), lo_chop(p%rcart(:,i),lo_sqtol), lo_chop(p%v(:,i),lo_sqtol) enddo end subroutine + + !> Writes a structure to QE output file or stdout + subroutine writetofile_qe(p,filename,write_velocities) + !> crystal structure + class(lo_crystalstructure), intent(in) :: p + !> filename + character(len=*), intent(in) :: filename + !> should be written. Default false. + logical, intent(in), optional :: write_velocities + ! local + integer :: i + + real(flyt) :: latpar + character(len=1000) :: opf + integer :: u + + ! file or stdout + if ( filename .eq. 'stdout' ) then + u=0 + else + u=open_file('out',filename) + endif + write(u,*) " nat = ", tochar(p%na) + write(u,*) " ntyp = ", tochar(p%nelements) + write(u,*) "/" + write(u,*) "CELL_PARAMETERS bohr" + write(u,*) lo_choplarge(p%latticevectors(:,1),lo_sqtol) + write(u,*) lo_choplarge(p%latticevectors(:,2),lo_sqtol) + write(u,*) lo_choplarge(p%latticevectors(:,3),lo_sqtol) + write(u,*) " " + write(u,*) "ATOMIC_POSITIONS bohr" + opf="(1X,a,2X,3F18.9))" + do i=1, p%na + write(u,opf) trim(p%atomic_symbol(p%species(i))), lo_chop(p%rcart(:,i),lo_sqtol) + enddo + if ( filename .ne. 'stdout' ) close(u) + end subroutine + + !> Writes a structure in PARSEC format + subroutine writetofile_parsec(p,filename,write_velocities) + !> crystal structure + class(lo_crystalstructure), intent(in) :: p + !> filename + character(len=*), intent(in) :: filename + !> should velocities be written. Default false. + logical, intent(in), optional :: write_velocities + ! local + integer :: i + ! + real(flyt) :: latpar + logical :: vel + integer :: u + + ! I don't have a velocity input format in parsec yet, so later problem + if ( present(write_velocities) ) then + vel=write_velocities + else + vel=.false. + endif + + ! file or stdout + if ( filename .eq. 'stdout' ) then + u=0 + else + u=open_file('out',filename) + endif + + write(u,*) "# comment line" + write(u,*) "Bulk" + if ( p%info%unitcell_lattice_parameter .gt. 0.0_flyt ) then + latpar=p%info%unitcell_lattice_parameter + write(u,"(2X,F21.13)") latpar*lo_bohr_to_A + write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,1)/latpar,1E-12_flyt) + write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,2)/latpar,1E-12_flyt) + write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,3)/latpar,1E-12_flyt) + else + latpar=1.0_flyt + write(u,"(2X,F20.12)") latpar + write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,1)/latpar,1E-12_flyt) + write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,2)/latpar,1E-12_flyt) + write(u,"(1X,3(1X,F22.14))") lo_chop(p%latticevectors(:,3)/latpar,1E-12_flyt) + endif + do i=1,p%na + write(u,"(1X,3(1X,E19.12),1X,A)") p%r(:,i),trim(p%atomic_symbol( p%species(i) )) + enddo + + if ( filename .ne. 'stdout' ) close(u) + end subroutine + end subroutine !> Read isotope distribution from file @@ -935,7 +1028,6 @@ module subroutine write_structure_to_hdf5(p,filename,input_id) !> in case I want to write it as a part of another file integer(HID_T), intent(in), optional :: input_id - integer :: i character(len=2000) :: atomnames type(lo_hdf5_helper) :: h5 diff --git a/src/libolle/type_forceconstant_secondorder.f90 b/src/libolle/type_forceconstant_secondorder.f90 index 76aca7b5..daa4de64 100644 --- a/src/libolle/type_forceconstant_secondorder.f90 +++ b/src/libolle/type_forceconstant_secondorder.f90 @@ -157,6 +157,8 @@ module type_forceconstant_secondorder procedure :: fake_forceconstant !> dump the forceconstant in Abinit anaddb format procedure :: write_to_anaddb + !> dump forceconstants into QE format + procedure :: write_to_qe !> Set parameters for Ewald summation and enforce hermiticity of Born effective charges procedure :: set_ewald_and_enforce_borncharge_hermiticity !> Return a dynamical matrix @@ -171,6 +173,8 @@ module type_forceconstant_secondorder procedure :: supercell_longrange_dynamical_matrix_at_gamma !> Get the commensurate modes for a supercell procedure :: commensurate_modes + !> + procedure :: write_dynmat_to_qe !> Size in memory, in bytes procedure :: size_in_mem => fc2_size_in_mem end type @@ -286,6 +290,23 @@ module subroutine write_to_anaddb(fc, uc, qgrid, mw, mem) type(lo_mpi_helper), intent(inout) :: mw type(lo_mem_helper), intent(inout) :: mem end subroutine + module subroutine write_to_qe(fc,uc,mw,mem) + !> forceconstant + class(lo_forceconstant_secondorder), intent(inout) :: fc + !> unitcell + type(lo_crystalstructure), intent(inout) :: uc + type(lo_mpi_helper), intent(inout) :: mw + type(lo_mem_helper), intent(inout) :: mem + end subroutine + module subroutine write_dynmat_to_qe(fc,uc,qgrid,mw,mem) + !> forceconstant + class(lo_forceconstant_secondorder), intent(inout) :: fc + !> unitcell + type(lo_crystalstructure), intent(inout) :: uc + integer, dimension(3),intent(in) :: qgrid + type(lo_mpi_helper), intent(inout) :: mw + type(lo_mem_helper), intent(inout) :: mem + end subroutine end interface contains diff --git a/src/libolle/type_forceconstant_secondorder_io.f90 b/src/libolle/type_forceconstant_secondorder_io.f90 index 84dda64f..76784712 100644 --- a/src/libolle/type_forceconstant_secondorder_io.f90 +++ b/src/libolle/type_forceconstant_secondorder_io.f90 @@ -1,7 +1,8 @@ submodule(type_forceconstant_secondorder) type_forceconstant_secondorder_io use konstanter, only: lo_forceconstant_2nd_HartreeBohr_to_eVA, lo_forceconstant_2nd_eVA_to_HartreeBohr, lo_emu_to_amu, & - lo_exitcode_io, lo_bohr_to_A + lo_exitcode_io, lo_bohr_to_A, lo_emu_to_amu use gottochblandat, only: open_file +use geometryfunctions, only: lo_inscribed_sphere_in_box use type_qpointmesh, only: lo_qpoint_mesh,lo_generate_qmesh use type_voronoi_distancetable, only: lo_voronoi_distancetable use mpi_wrappers, only: lo_mpi_helper @@ -266,77 +267,77 @@ module subroutine write_to_anaddb(fc,uc,qgrid,mw,mem) ! Abinit wants dynamical matrices, I have forceconstants. I need to generate a q-mesh ! tight enough to ensure that there is no information loss. Will probably pad it a little ! to be on the safe side. -! findqmesh: block -! type(lo_voronoi_distancetable) :: dt -! real(r8), dimension(3) :: v0 -! real(r8) :: cutoff -! integer, dimension(:,:), allocatable :: supercells -! integer, dimension(3) :: dims -! integer :: i,j,k,l,nq -! -! ! I always screw up the cutoff, do it again here to be on the safe side -! cutoff=0.0_r8 -! do i=1,fc%na -! do j=1,fc%atom(i)%n -! cutoff=max(cutoff,norm2(fc%atom(i)%pair(j)%r)) -! enddo -! enddo -! ! Get a bunch of reasonable supercells -!! call generate_possible_supercells(uc,cutoff,supercells) -!! ! Find one large enough to fit the forceconstant -!! ssl1: do i=1,size(supercells,2) -!! call dt%generate_wzcutoff( uc%r,uc%latticevectors,supercells(:,i) ) -!! if ( does_forceconstant_fit(fc,dt) ) then -!! dims=supercells(:,i) -!! exit ssl1 -!! endif -!! enddo ssl1 -! dims=(/6,6,6/) -! ! add a little safety margin to the dimensions, not like it actually costs anything -! dims=dims+2 -! -! ! build the list of qvectors -! nq=product(dims) -! allocate(qvectors(3,nq)) -! l=0 -! do i=1,dims(1) -! do j=1,dims(2) -! do k=1,dims(3) -! l=l+1 -! v0(1)=(i-1.0_r8)/(1.0_r8*dims(1)) -! v0(2)=(j-1.0_r8)/(1.0_r8*dims(2)) -! v0(3)=(k-1.0_r8)/(1.0_r8*dims(3)) -! qvectors(:,l)=uc%fractional_to_cartesian(v0,reciprocal=.true.) -! enddo -! enddo -! enddo -! end block findqmesh + ! findqmesh: block + ! type(lo_voronoi_distancetable) :: dt + ! real(r8), dimension(3) :: v0 + ! real(r8) :: cutoff + ! integer, dimension(:,:), allocatable :: supercells + ! integer, dimension(3) :: dims + ! integer :: i,j,k,l,nq + ! + ! ! I always screw up the cutoff, do it again here to be on the safe side + ! cutoff=0.0_r8 + ! do i=1,fc%na + ! do j=1,fc%atom(i)%n + ! cutoff=max(cutoff,norm2(fc%atom(i)%pair(j)%r)) + ! enddo + ! enddo + ! ! Get a bunch of reasonable supercells + !! call generate_possible_supercells(uc,cutoff,supercells) + !! ! Find one large enough to fit the forceconstant + !! ssl1: do i=1,size(supercells,2) + !! call dt%generate_wzcutoff( uc%r,uc%latticevectors,supercells(:,i) ) + !! if ( does_forceconstant_fit(fc,dt) ) then + !! dims=supercells(:,i) + !! exit ssl1 + !! endif + !! enddo ssl1 + ! dims=(/6,6,6/) + ! ! add a little safety margin to the dimensions, not like it actually costs anything + ! dims=dims+2 + ! + ! ! build the list of qvectors + ! nq=product(dims) + ! allocate(qvectors(3,nq)) + ! l=0 + ! do i=1,dims(1) + ! do j=1,dims(2) + ! do k=1,dims(3) + ! l=l+1 + ! v0(1)=(i-1.0_r8)/(1.0_r8*dims(1)) + ! v0(2)=(j-1.0_r8)/(1.0_r8*dims(2)) + ! v0(3)=(k-1.0_r8)/(1.0_r8*dims(3)) + ! qvectors(:,l)=uc%fractional_to_cartesian(v0,reciprocal=.true.) + ! enddo + ! enddo + ! enddo + ! end block findqmesh ! Now I have a decent list of q-vectors. For each of those, I need a dynamical matrix. getdynmat: block integer :: i,nb complex(r8), dimension(:,:,:), allocatable :: Dq - real(r8) :: abiqpts(3,12) + !real(r8) :: abiqpts(3,12) nb=uc%na*3 call lo_generate_qmesh(qp,uc,qgrid,'fft',timereversal=.true.,headrankonly=.false.,& mw=mw,mem=mem,verbosity=1) nq = qp%n_irr_point -!nq = 12 -!abiqpts =reshape ((/ & -! 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, & -! 2.50000000E-01, 0.00000000E+00, 0.00000000E+00, & -! 5.00000000E-01, 0.00000000E+00, 0.00000000E+00, & -! 2.50000000E-01, 2.50000000E-01, 0.00000000E+00, & -! 0.00000000E+00, 0.00000000E+00, 2.50000000E-01, & -! 2.50000000E-01, 0.00000000E+00, 2.50000000E-01, & -! 5.00000000E-01, 0.00000000E+00, 2.50000000E-01, & -! 2.50000000E-01, 2.50000000E-01, 2.50000000E-01, & -! 0.00000000E+00, 0.00000000E+00, 5.00000000E-01, & -! 2.50000000E-01, 0.00000000E+00, 5.00000000E-01, & -! 5.00000000E-01, 0.00000000E+00, 5.00000000E-01, & -! 2.50000000E-01, 2.50000000E-01, 5.00000000E-01 /), (/3,12/) ) + !nq = 12 + !abiqpts =reshape ((/ & + ! 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, & + ! 2.50000000E-01, 0.00000000E+00, 0.00000000E+00, & + ! 5.00000000E-01, 0.00000000E+00, 0.00000000E+00, & + ! 2.50000000E-01, 2.50000000E-01, 0.00000000E+00, & + ! 0.00000000E+00, 0.00000000E+00, 2.50000000E-01, & + ! 2.50000000E-01, 0.00000000E+00, 2.50000000E-01, & + ! 5.00000000E-01, 0.00000000E+00, 2.50000000E-01, & + ! 2.50000000E-01, 2.50000000E-01, 2.50000000E-01, & + ! 0.00000000E+00, 0.00000000E+00, 5.00000000E-01, & + ! 2.50000000E-01, 0.00000000E+00, 5.00000000E-01, & + ! 5.00000000E-01, 0.00000000E+00, 5.00000000E-01, & + ! 2.50000000E-01, 2.50000000E-01, 5.00000000E-01 /), (/3,12/) ) allocate(dynmat(nb,nb,nq)) allocate(Dq(nb,nb,3)) @@ -344,9 +345,9 @@ module subroutine write_to_anaddb(fc,uc,qgrid,mw,mem) do i=1,nq qvectors(:,i) = uc%cartesian_to_fractional(qp%ip(i)%r,reciprocal=.true.) -!qvectors(:,i) = abiqpts(:,i) -!qp%ip(i)%r = uc%fractional_to_cartesian(abiqpts(:,i), reciprocal=.true.) -!print *, 'iq qpt ', i, qvectors(:,i) + !qvectors(:,i) = abiqpts(:,i) + !qp%ip(i)%r = uc%fractional_to_cartesian(abiqpts(:,i), reciprocal=.true.) + !print *, 'iq qpt ', i, qvectors(:,i) call fc%dynamicalmatrix(uc,qp%ip(i),dynmat(:,:,i),mem,Dq,qdirection=(/1.0d0,0.0d0,0.0d0/),skipnonanalytical=.true.) enddo deallocate(Dq) @@ -356,7 +357,7 @@ module subroutine write_to_anaddb(fc,uc,qgrid,mw,mem) printstuff: block real(r8), dimension(:,:), allocatable :: tnons real(r8), dimension(:,:), allocatable :: spinat - real(r8), dimension(:, :), allocatable :: rotmat + !real(r8), dimension(:, :), allocatable :: rotmat real(r8), dimension(:), allocatable :: amu,znucl ! dumna,dummass integer, dimension(:,:,:), allocatable :: symrel integer, dimension(:), allocatable :: symafm,typat @@ -369,18 +370,17 @@ module subroutine write_to_anaddb(fc,uc,qgrid,mw,mem) &9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 11,& &12, 13, 14, 15, 16, 17, 18, 19, 12, 13, 14, 15, 16, 7, 8,& &9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 11] -! &'H ', 'He',& -! &'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', 'Ne',& -! &'Na', 'Mg', 'Al', 'Si', 'P ', 'S ', 'Cl', 'Ar',& -! &'K ', 'Ca', 'Sc', 'Ti', 'V ', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr',& -! &'Rb', 'Sr', 'Y ', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I ', 'Xe',& -! &'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', & -! & 'Hf', 'Ta', 'W ', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn',& -! &'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U ', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr'] + ! &'H ', 'He',& + ! &'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', 'Ne',& + ! &'Na', 'Mg', 'Al', 'Si', 'P ', 'S ', 'Cl', 'Ar',& + ! &'K ', 'Ca', 'Sc', 'Ti', 'V ', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr',& + ! &'Rb', 'Sr', 'Y ', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I ', 'Xe',& + ! &'Cs', 'Ba', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', & + ! & 'Hf', 'Ta', 'W ', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn',& + ! &'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U ', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr'] integer, dimension(18) :: ngfft - integer :: i,j,a1,a2,ii,jj,x,y,nb - integer :: u, nq + integer :: i,j,nb integer :: ntypat character(len=500) :: string character(len=500) :: filnam @@ -476,13 +476,13 @@ module subroutine write_to_anaddb(fc,uc,qgrid,mw,mem) end block printstuff - printdm : block + printdm : block complex(r8), dimension(:,:), allocatable :: dmt real(r8), dimension(:,:), allocatable :: rotmat real(r8), dimension(3,3) :: bec real(r8), dimension(3,3) :: eps real(r8), dimension(3,3) :: m1,m0 - integer :: i,j,a1,a2,ii,jj,x,y,nb + integer :: i,a1,a2,ii,jj,x,y,nb integer :: npert m1=0.0d0 @@ -583,6 +583,384 @@ module subroutine write_to_anaddb(fc,uc,qgrid,mw,mem) deallocate(qvectors) end subroutine +module subroutine write_to_qe(fc,uc,mw,mem) + !> forceconstant + class(lo_forceconstant_secondorder), intent(inout) :: fc + !> unitcell + type(lo_crystalstructure), intent(inout) :: uc + type(lo_mpi_helper), intent(inout) :: mw + type(lo_mem_helper), intent(inout) :: mem + + !type(lo_forceconstant_secondorder) :: fcss + real(r8), dimension(:,:,:,:,:,:,:), allocatable :: rm1 + integer, dimension(3) :: supercelldim + + init: block + real(r8), parameter :: safety_margin_for_cutoff=0.1_r8 + type(lo_crystalstructure) :: ss + + real(r8), dimension(3,3) :: m0 + real(r8), dimension(3) :: v0 !,v1,v2 + real(r8) :: f0 + integer, dimension(3) :: gi,gj + integer :: i,j,k,a1,a2,ii,jj + + ! So ... first step would be to find a (diagonal) supercell that safely bounds + ! the TDEP forceconstant with some margin. + + + supercelldim=1 + do + do i = 1, 3 + m0(:, i) = uc%latticevectors(:,i)*supercelldim(i) + end do + f0 = lo_inscribed_sphere_in_box(m0) + if (f0 .gt. fc%cutoff+safety_margin_for_cutoff ) exit + supercelldim = increment_dimensions(supercelldim, uc%latticevectors) + end do + supercelldim=supercelldim+2 + + ! build the supercell + call uc%build_supercell(ss,dimensions=supercelldim) + ! establish mapping between unit and supercell + call ss%classify('supercell', uc) + ! remap the forceconstant so that it is defined for the supercell + !call fc%remap(uc, ss, fcss) + + ! I might want to have this in (xyz,xyz,n_a,sx,sy,sz) instead. + ! So, let's give this a shot. Might have to think a little. + allocate(rm1(3,3,uc%na,uc%na,supercelldim(1),supercelldim(2),supercelldim(3))) + rm1=0.0_r8 + + k=0 ! counter for matches + do i=1,ss%na + gi=ss%info%cellindex(:,i) + if ( sum(abs(gi-[1,1,1])) .ne. 0 ) cycle + a1=ss%info%index_in_unitcell(i) + do j=1,ss%na + a2=ss%info%index_in_unitcell(j) + gj=ss%info%cellindex(:,j) + + ! vector between atoms, fractional supercell coordinates + v0=ss%r(:,j)-ss%r(:,i) + !v0=ss%r(:,i)-ss%r(:,j) + ! adjust with pbc? Think so + v0=lo_clean_fractional_coordinates(v0+0.5_r8)-0.5_r8 + ! to Cartesian + v0=matmul(ss%latticevectors,v0) + + ! This should be matched to a TDEP ifc! + + do ii=1,fc%atom(a1)%n + if ( fc%atom(a1)%pair(ii)%i2 .ne. a2 ) cycle + if ( norm2(v0 -fc%atom(a1)%pair(ii)%r) .lt. 1E-5_r8 ) then + k=k+1 + ! we have a match! This should be enough I think. + !rm1(:,:,a1,a2,gj(1),gj(2),gj(3)) = transpose( fc%atom(a1)%pair(ii)%m ) + rm1(:,:,a2,a1,gj(1),gj(2),gj(3)) = transpose( fc%atom(a1)%pair(ii)%m ) + endif + enddo + enddo + enddo + + ! Sanity check that all pairs got matched. Or at least the correct numbers + jj=0 + do a1=1,fc%na + jj=jj+fc%atom(a1)%n + enddo + if ( jj .ne. k ) call lo_stop_gracefully(['Did not match all pairs'],lo_exitcode_symmetry,__FILE__,__LINE__) + + + end block init + + if ( mw%talk ) then + dumpstuff: block + character(len=1000) :: str0,str1 + real(r8) :: f0 + integer :: u,i,j,a1,a2,ii,jj,kk + + u = open_file('out','outfile.qe_fc.xml') + ! Let's start with the header + write(u,"(A)") '' + write(u,"(A)") '' + write(u,"(A)") ' ' + write(u,"(A)") ' '//tochar(uc%nelements)//'' + write(u,"(A)") ' '//tochar(uc%na)//'' + write(u,"(A)") ' 0' + write(u,"(A)") ' 1' + write(u,"(A)") ' ' + write(u,"(A)") ' 1.0E+00 0.000000000000000E+00 0.000000000000000E+00' + write(u,"(A)") ' 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00' + write(u,"(A)") ' ' + write(u,"(A)") ' ' + write(u,*) uc%latticevectors(:,1) + write(u,*) uc%latticevectors(:,2) + write(u,*) uc%latticevectors(:,3) + write(u,"(A)") ' ' + write(u,"(A)") ' ' + write(u,*) uc%reciprocal_latticevectors(:,1) + write(u,*) uc%reciprocal_latticevectors(:,2) + write(u,*) uc%reciprocal_latticevectors(:,3) + write(u,"(A)") ' ' + write(u,"(A)") ' '//tochar(uc%volume)//'' + do i=1,uc%nelements + write(u,"(A)") ' '//trim(uc%atomic_symbol(i))//'' + f0=-1.0_r8 + do j=1,uc%na + if ( uc%species(j) .eq. i ) then + f0=uc%isotope(j)%mean_mass*lo_emu_to_amu + exit + endif + enddo + + write(u,"(A)") ' '//trim( tochar(f0,13) )//'' + enddo + do i=1,uc%na + j=uc%species(i) + str0=' ' + write(u,"(A)") trim(str0) + enddo + write(u,"(A)") ' 6' + write(u,"(A)") ' ' + + ! post header is dielectric stuff + if ( fc%polar ) then + write(u,"(A)") ' ' + write(u,"(A)") ' ' + do i=1,3 + write(u,*) fc%loto%eps(i,:) + enddo + write(u,"(A)") ' ' + write(u,"(A)") ' ' + do i=1,fc%na + write(u,"(A)") ' ' + do j=1,3 + write(u,*) fc%loto%born_effective_charges(j,:,i) !*0.5_r8 + enddo + write(u,"(A)") ' ' + enddo + write(u,"(A)") ' ' + write(u,"(A)") ' ' + else + write(u,"(A)") ' ' + write(u,"(A)") ' ' + endif + + ! Now we start with the actual force constants! In some units and some convention. + write(u,"(A)") ' ' + + write(u,"(A)") ' ' + write(u,*) supercelldim + write(u,"(A)") ' ' + if ( fc%polar ) then + f0=fc%ew%lambda + f0=(f0/lo_twopi)**2 + ! Maybe it should be 1/, sqrt, or sth + write(u,"(A)") ' '//tochar(f0,12)//'' + else + write(u,"(A)") ' 1E-300' + endif + ! Now for the actual IFCs. + + do a1=1,uc%na + do a2=1,uc%na + do kk=1,supercelldim(3) + do jj=1,supercelldim(2) + do ii=1,supercelldim(1) + + str0="s_s1_m1_m2_m3."//tochar(a1)//"."//tochar(a2)//"."//tochar(ii)//"."//tochar(jj)//"."//tochar(kk) + write(u,"(A)") ' <'//trim(str0)//'>' + write(u,"(A)") ' ' + do i=1,3 + write(u,*) rm1(:,i,a1,a2,ii,jj,kk)*2.0_r8 ! Factor 2 because Ha->Ry + enddo + write(u,"(A)") ' ' + write(u,"(A)") ' ' + enddo + enddo + enddo + enddo + enddo + write(u,"(A)") ' ' + write(u,"(A)") '' + write(u,"(A)") '' + close(u) + end block dumpstuff + endif + +end subroutine + +module subroutine write_dynmat_to_qe(fc,uc,qgrid,mw,mem) + !> forceconstant + class(lo_forceconstant_secondorder), intent(inout) :: fc + !> unitcell + type(lo_crystalstructure), intent(inout) :: uc + !> q-grid density + integer, dimension(3),intent(in) :: qgrid + !> MPI helper + type(lo_mpi_helper), intent(inout) :: mw + !> memory tracker + type(lo_mem_helper), intent(inout) :: mem + + class(lo_qpoint_mesh), allocatable :: qp + + ! Generate the q-grid + call lo_generate_qmesh(qp,uc,qgrid,'fft',timereversal=.true.,headrankonly=.false.,mw=mw,mem=mem,verbosity=1) + + ! Write the header main xml file thingy? + writesummary: block + real(r8), dimension(3) :: v0 + integer :: u,i + + u = open_file('out','outfile.qe_dyn0') + write(u,*) qgrid + write(u,*) qp%n_full_point + do i=1,qp%n_full_point + v0=qp%ap(i)%r + write(u,*) v0 + enddo + close(u) + end block writesummary + + dumpdynmat: block + complex(r8), dimension(:,:), allocatable :: dynamical_matrix + complex(r8), dimension(3,3) :: cm0 + real(r8) :: f0 + integer :: iq,u,i,j,a1,a2,ii,jq + character(len=1000) :: str0,str1 + + allocate(dynamical_matrix(uc%na*3,uc%na*3)) + dynamical_matrix=0.0_r8 + + do iq=1,qp%n_full_point + u = open_file('out','outfile.qe_dyn'//tochar(iq)//'.xml') + write(u,"(A)") '' + write(u,"(A)") '' + write(u,"(A)") ' ' + write(u,"(A)") ' '//tochar(uc%nelements)//'' + write(u,"(A)") ' '//tochar(uc%na)//'' + write(u,"(A)") ' 0' + write(u,"(A)") ' 1' + write(u,"(A)") ' ' + write(u,"(A)") ' 1.0E+00 0.000000000000000E+00 0.000000000000000E+00' + write(u,"(A)") ' 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00' + write(u,"(A)") ' ' + write(u,"(A)") ' ' + write(u,*) uc%latticevectors(:,1) + write(u,*) uc%latticevectors(:,2) + write(u,*) uc%latticevectors(:,3) + write(u,"(A)") ' ' + write(u,"(A)") ' ' + write(u,*) uc%reciprocal_latticevectors(:,1) + write(u,*) uc%reciprocal_latticevectors(:,2) + write(u,*) uc%reciprocal_latticevectors(:,3) + write(u,"(A)") ' ' + write(u,"(A)") ' '//tochar(uc%volume)//'' + do i=1,uc%nelements + write(u,"(A)") ' '//trim(uc%atomic_symbol(i))//'' + f0=-1.0_r8 + do j=1,uc%na + if ( uc%species(j) .eq. i ) then + f0=uc%isotope(j)%mean_mass*lo_emu_to_amu + exit + endif + enddo + + write(u,"(A)") ' '//trim( tochar(f0,13) )//'' + enddo + do i=1,uc%na + j=uc%species(i) + str0=' ' + write(u,"(A)") trim(str0) + enddo + write(u,"(A)") ' 1' + write(u,"(A)") ' ' + if ( fc%polar ) then + write(u,"(A)") ' ' + write(u,"(A)") ' ' + do i=1,3 + write(u,*) fc%loto%eps(i,:) + enddo + write(u,"(A)") ' ' + write(u,"(A)") ' ' + do i=1,fc%na + write(u,"(A)") ' ' + do j=1,3 + write(u,*) fc%loto%born_effective_charges(j,:,i) !*0.5_r8 + enddo + write(u,"(A)") ' ' + enddo + write(u,"(A)") ' ' + write(u,"(A)") ' ' + else + write(u,"(A)") ' ' + write(u,"(A)") ' ' + endif + + ! And now to write the actual q-point information + write(u,"(A)") '' + + write(u,"(A)") '' + write(u,*) qp%ap(iq)%r + write(u,"(A)") '' + + ! Ok so this is a reasonable place to actually calculate the dynamical matrix. + call fc%dynamicalmatrix(uc,qp%ap(iq),dynamical_matrix,mem,skipnonanalytical=.true.) + + do a2=1,uc%na + do a1=1,uc%na + cm0=dynamical_matrix( (a1-1)*3+1:a1*3, (a2-1)*3+1:a2*3 ) + cm0=cm0/(uc%invsqrtmass(a1)*uc%invsqrtmass(a2)) + cm0=cm0*2.0_r8 !*0.5_r8 + write(u,"(A)") '' + do j=1,3 + do i=1,3 + write(u,*) real(cm0(i,j),r8),-aimag(cm0(i,j)) + enddo + enddo + write(u,"(A)") '' + enddo + enddo + + write(u,"(A)") '' + write(u,"(A)") '' + do i=1,uc%na*3 + write(u,"(A)") '' + write(u,*) 1.0_r8,2.0_r8 + write(u,"(A)") '' + write(u,"(A)") '' + do j=1,uc%na*3 + write(u,*) 3.0_r8,4.0_r8 + enddo + write(u,"(A)") '' + enddo + write(u,"(A)") '' + + write(u,"(A)") '' + write(u,"(A)") '' + close(u) + enddo + end block dumpdynmat + + dumpq2r: block + integer :: u + + u = open_file('out','outfile.qe_q2r') + write(u,"(A)") '&input' + write(u,"(A)") "zasr='simple', fildyn='outfile.qe_dyn.xml', flfrc='qe_fc'" + write(u,"(A)") "/" + close(u) + end block dumpq2r + + write(*,*) 'DONE HERE FOR NOW' + stop + +end subroutine + !> mini-version of my dynamical matrix calculator, repeated here to avoid circular dependencies subroutine simple_dynmat(fc, uc, qv, dynmat) !> forceconstant @@ -933,4 +1311,49 @@ subroutine ddb_io_out(dscrpt, filnam, matom, mband,& end do end subroutine ddb_io_out +function increment_dimensions(dimin, box) result(dimut) + integer, dimension(3), intent(in) :: dimin + real(r8), dimension(3, 3), intent(in) :: box + integer, dimension(3) :: dimut + ! + real(r8), dimension(3, 3) :: m0 + integer, dimension(3) :: di + integer :: i, j, k + real(r8) :: f0, f1, ff0 + + ! try with the sphere thing. First get a baseline + do j = 1, 3 + m0(:, j) = box(:, j)*(2*dimin(j) + 1) + end do + ff0 = lo_inscribed_sphere_in_box(m0) + ! Increment the dimension that gives the biggest increase in radii + f0 = 0.0_r8 + dimut = 0 + do i = 1, 3 + di = dimin + di(i) = di(i) + 1 + do j = 1, 3 + m0(:, j) = box(:, j)*(2*di(j) + 1) + end do + f1 = lo_inscribed_sphere_in_box(m0) + if (f1 .gt. f0 .and. abs(f1 - ff0) .gt. lo_tol) then + dimut = di + f0 = f1 + end if + end do + + ! if nothing helped, increment the lowest number + if (dimut(1) .eq. 0) then + j = lo_hugeint + do i = 1, 3 + if (di(i) .lt. j) then + j = di(i) + k = i + end if + end do + dimut = dimin + dimut(k) = dimut(k) + 1 + end if +end function + end submodule