From a23801ff0da9b805a4abc63ac81e842371391da9 Mon Sep 17 00:00:00 2001 From: David Morse Date: Fri, 16 Aug 2019 21:44:53 -0500 Subject: [PATCH] Adding code to output groups in pscfpp format. --- src/group_mod.f90 | 231 +++++++++++++++++++++++++++++++++++++++------- src/pscf.fp.f90 | 15 +++ 2 files changed, 212 insertions(+), 34 deletions(-) diff --git a/src/group_mod.f90 b/src/group_mod.f90 index 2fb85ae..fe4dbfc 100644 --- a/src/group_mod.f90 +++ b/src/group_mod.f90 @@ -101,7 +101,6 @@ module group_mod real(long) :: v(3) ! translation vector end type symmetry_type !*** - !------------------------------------------------------------------ !****t group_mod/group_type @@ -318,66 +317,230 @@ end subroutine output_matrix !=================================================================== + !------------------------------------------------------------------- !****p group_mod/output_symmetry ! SUBROUTINE - ! output_symmetry(s,iunit) + ! output_symmetry(s, iunit, mode_arg) ! PURPOSE ! Write symmetry s to file iunit ! SOURCE !------------------------------------------------------------------- - subroutine output_symmetry(s,iunit) - type(symmetry_type), intent(in):: s - integer :: iunit + subroutine output_symmetry(s, iunit, mode_arg) + type(symmetry_type), intent(in) :: s + integer :: iunit + integer, intent(in), optional :: mode_arg !*** - integer :: i, j + integer :: i, j, mode + + ! Set output format mode + if (present(mode_arg)) then + mode = mode_arg + else + mode = 1 + end if + + if (mode == 1) then + + ! Output in native (fortran pscf) format + do i=1, dim + if (dim == 3) then + write(iunit,FMT='(3F8.3,5X,F8.3)') & + (s%m(i,j),j=1, dim), s%v(i) + else if (dim == 2) then + write(iunit,FMT='(2F8.3,5X,F8.3)') & + (s%m(i,j),j=1, dim), s%v(i) + else if (dim == 1) then + write(iunit,FMT='(F8.3,5X,F8.3)') & + (s%m(i,j),j=1, dim), s%v(i) + else + write(iunit,*) 'Invalid dim in output_symmetry' + endif + enddo + write(iunit,*) + + else if (mode == 2) then + + ! Output in pscfpp (c++ code) format + call output_symmetry_c(s, iunit); + + end if + + end subroutine output_symmetry + !=================================================================== + + !------------------------------------------------------------------- + !****p group_mod/output_symmetry_c + ! SUBROUTINE + ! output_symmetry_c(s,iunit) + ! PURPOSE + ! Write symmetry s to file iunit in format of C++ pscfpp code. + ! SOURCE + !------------------------------------------------------------------- + subroutine output_symmetry_c(s, iunit) + type(symmetry_type), intent(in) :: s + integer :: iunit + !*** + real(long) :: r, q + integer :: i, j, k, d, n + integer :: rot(3,3) ! point group matrix + integer :: den(3) ! denominator of translation fraction + integer :: num(3) ! numerator of translation fraction + integer :: trial(4) ! allowed values of denominator + logical :: found + + ! Convert elements of matrix m to nearest integers do i=1, dim - if (dim == 3) then - write(iunit,FMT='(3F8.3,5X,F8.3)') & - (s%m(i,j),j=1, dim),s%v(i) - else if (dim == 2) then - write(iunit,FMT='(2F8.3,5X,F8.3)') & - (s%m(i,j),j=1, dim),s%v(i) - else if (dim == 1) then - write(iunit,FMT='(F8.3,5X,F8.3)') & - (s%m(i,j),j=1, dim),s%v(i) - else - write(iunit,*) 'invalid dim in output_symmetry' - endif + do j=1, dim + r = anint(s%m(i, j)) + if (dabs(r - s%m(i,j)) > epsilon) then + write(iunit,*) 'Rotation matrix element is non-integer' + return + else + rot(i, j) = nint(r) + endif + enddo + enddo + + ! Create allowed values of denominator + trial(1) = 2 + trial(2) = 3 + trial(3) = 4 + trial(4) = 6 + + ! Convert elements of vector s%v to nearest fractions + do i=1, dim + found = .false. + r = anint(s%v(i)) + if (abs(r - s%v(i)) < 1.0E-3) then + found = .true. + num(i) = 0 + den(i) = 1 + else + do j=1, 4 + k = trial(j) + q = k*s%v(i) + r = anint(q) + if (abs(r - q) < 1.0E-3) then + found = .true. + den(i) = k + num(i) = nint(q) + q = real(num(i))/real(den(i)) + if (abs(q - s%v(i)) > 1.0E-2) then + write(iunit,*) 'Incorrect conversion to fraction' + write(iunit,*) 'Translation element = ', s%v(i) + write(iunit,*) 'Numerator = ', num(i) + write(iunit,*) 'Denominator = ', den(i) + return + endif + exit + endif + enddo + if (den(i) < 2) then + write(iunit,*) 'Invalid denominator' + write(iunit,*) 'Translation element = ', s%v(i) + write(iunit,*) 'Numerator = ', num(i) + write(iunit,*) 'Denominator = ', den(i) + return + endif + endif + + if (.not.found) then + write(iunit,*) 'Translation element = ', s%v(i) + write(iunit,*) 'Translation element is not small fraction' + return + endif + + enddo + + ! Write rotation matrix + do i=1, dim + if (dim == 3) then + write(iunit,FMT='(3I6)') (rot(i,j), j=1, dim) + else if (dim == 2) then + write(iunit,FMT='(2I6)') (rot(i,j), j=1, dim) + else if (dim == 1) then + write(iunit,FMT='(I6)') (rot(i,j), j=1, dim) + else + write(iunit,*) 'Invalid dim in output_symmetry' + endif + enddo + + ! Write translation vector + do i=1, dim + if (num(i) == 0) then + write(iunit,FMT='(I6)', advance="no") 0 + else + write(iunit,FMT='(I4,"/",I1)', advance="no") num(i), den(i) + endif enddo write(iunit,*) - end subroutine output_symmetry + + write(iunit,*) + end subroutine output_symmetry_c !=================================================================== !------------------------------------------------------------------- !****p group_mod/output_group ! SUBROUTINE - ! output_group(g,iunit) + ! output_group(g,iunit,mode_arg) ! PURPOSE ! Write group g to file iunit ! SOURCE !------------------------------------------------------------------- - subroutine output_group(g,iunit) - type(group_type), intent(IN) :: g - integer, intent(IN) :: iunit + subroutine output_group(g,iunit,mode_arg) + type(group_type), intent(IN) :: g + integer, intent(IN) :: iunit + integer, intent(IN), optional :: mode_arg !*** type(version_type) :: version + integer :: mode integer :: i - version%major = 1 - version%minor = 0 - call output_version(version, iunit) - write(iunit,FMT="(A9,' = basis')") g%s(1)%basis - write(iunit,FMT="(i9,' = dimension')") dim - write(iunit,FMT="(i9,' = # of symmetries')") g%order - write(iunit,*) - do i=1, g%order - write(iunit,FMT='(i3)') i - call output_symmetry(g%s(i),iunit) - enddo + ! Check for optional mode argument + if (present(mode_arg)) then + mode = mode_arg + else + mode = 1 + end if + + ! Output preamble (if any) + if (mode == 1) then + + ! Write preamble + version%major = 1 + version%minor = 0 + call output_version(version, iunit) + write(iunit,FMT="(A9,' = basis')") g%s(1)%basis + write(iunit,FMT="(i9,' = dimension')") dim + write(iunit,FMT="(i9,' = # of symmetries')") g%order + write(iunit,*) + + ! Write symmetry elements + do i=1, g%order + write(iunit,FMT='(i3)') i + call output_symmetry(g%s(i),iunit, mode) + enddo + + else if (mode == 2) then + + ! Write number of symmetry elements + write(iunit,FMT="(i9)") g%order + write(iunit,*) + + ! Write symmetry elements + do i=1, g%order + call output_symmetry(g%s(i),iunit, mode) + enddo + + else + + write(iunit, *) "Incorrect mode format number = ", mode + + endif end subroutine output_group !=================================================================== diff --git a/src/pscf.fp.f90 b/src/pscf.fp.f90 index 987e89d..450f95d 100644 --- a/src/pscf.fp.f90 +++ b/src/pscf.fp.f90 @@ -735,6 +735,21 @@ program pscf call output_group(group,field_unit) close(field_unit) + case ('OUTPUT_GROUP_CPP') + + ! Check preconditions (Needs group created in BASIS block) + if (.not. basis_flag) then + write(6,*) "Error: Must read BASIS before OUTPUT_GROUP" + exit op_loop + end if + iterate_flag = .FALSE. + omega_flag = .FALSE. + + call input(output_filename,'output_filename') + open(unit=field_unit,file=trim(output_filename), status='replace') + call output_group(group,field_unit, 2) + close(field_unit) + case ('OUTPUT_WAVES') ! Check preconditions (Needs BASIS)