w90_readwrite_get_projections Subroutine
public subroutine w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann, input_proj, inv_lattice, lcount, spinors, bohr, stdout, error, comm)
Fills the projection data block
Arguments
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(settings_type), | intent(inout) | :: | settings | |||
| integer, | intent(inout) | :: | num_proj | |||
| type(atom_data_type), | intent(in) | :: | atom_data | |||
| integer, | intent(in) | :: | num_wann | |||
| type(proj_type), | intent(inout), | allocatable | :: | input_proj(:) | ||
| real(kind=dp), | intent(in) | :: | inv_lattice(3,3) | |||
| logical, | intent(in) | :: | lcount | |||
| logical, | intent(in) | :: | spinors | |||
| real(kind=dp), | intent(in) | :: | bohr | |||
| integer, | intent(in) | :: | stdout | |||
| type(w90_error_type), | intent(out), | allocatable | :: | error | ||
| type(w90_comm_type), | intent(in) | :: | comm |
Calls
Called by
Source Code
subroutine w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann, input_proj, & inv_lattice, lcount, spinors, bohr, stdout, error, comm) !================================================! ! !! Fills the projection data block ! !================================================! use w90_constants, only: eps6, eps2 use w90_utility, only: utility_cart_to_frac, utility_string_to_coord, utility_strip use w90_error, only: w90_error_type, set_error_alloc, set_error_input implicit none ! arguments type(atom_data_type), intent(in) :: atom_data type(proj_type), allocatable, intent(inout) :: input_proj(:) type(w90_error_type), allocatable, intent(out) :: error type(w90_comm_type), intent(in) :: comm type(settings_type), intent(inout) :: settings integer, intent(in) :: num_wann integer, intent(inout) :: num_proj integer, intent(in) :: stdout real(kind=dp), intent(in) :: bohr real(kind=dp), intent(in) :: inv_lattice(3, 3) logical, intent(in) :: lcount logical, intent(in) :: spinors ! local variables real(kind=dp) :: pos_frac(3) real(kind=dp) :: pos_cart(3) character(len=20) :: keyword integer :: in, ins, ine, loop, line_e, line_s, block_line_e, block_line_s, counter integer :: sites, species, line, pos1, pos2, pos3, m_tmp, l_tmp, mstate integer :: loop_l, loop_m, loop_sites, ierr, loop_s, spn_counter logical :: found_e, found_s character(len=maxlen) :: dummy, end_st, start_st character(len=maxlen) :: ctemp, ctemp2, ctemp3, ctemp4, ctemp5, m_string integer, parameter :: min_l = -5 integer, parameter :: max_l = 3 integer, parameter :: min_m = 1 integer, parameter :: max_m = 7 integer :: ang_states(min_m:max_m, min_l:max_l) ! default values for the optional part of the projection definitions real(kind=dp), parameter :: proj_z_def(3) = (/0.0_dp, 0.0_dp, 1.0_dp/) real(kind=dp), parameter :: proj_x_def(3) = (/1.0_dp, 0.0_dp, 0.0_dp/) real(kind=dp), parameter :: proj_s_qaxis_def(3) = (/0.0_dp, 0.0_dp, 1.0_dp/) real(kind=dp), parameter :: proj_zona_def = 1.0_dp integer, parameter :: proj_radial_def = 1 real(kind=dp) :: proj_z_tmp(3) real(kind=dp) :: proj_x_tmp(3) real(kind=dp) :: proj_s_qaxis_tmp(3) real(kind=dp) :: proj_zona_tmp integer :: proj_radial_tmp logical :: lconvert, lrandom, proj_u_tmp, proj_d_tmp, found_f logical :: lpartrandom real(kind=dp) :: xnorm, znorm, cosphi, sinphi, xnorm_new, cosphi_new keyword = "projections" found_s = .false. found_e = .false. start_st = 'begin '//trim(keyword) end_st = 'end '//trim(keyword) ierr = 0 if (allocated(input_proj)) return ! projectors have already been read, return if (.not. lcount) then allocate (input_proj(num_proj), stat=ierr) if (ierr /= 0) then call set_error_alloc(error, 'Error allocating input_proj in w90_readwrite_get_projections', comm) return end if end if lconvert = .false. lrandom = .false. lpartrandom = .false. if (allocated(settings%in_data)) then ! we are reading from the input file do loop = 1, settings%num_lines ins = index(settings%in_data(loop), trim(keyword)) if (ins == 0) cycle in = index(settings%in_data(loop), 'begin') if (in == 0 .or. in > 1) cycle line_s = loop if (found_s) then call set_error_input(error, 'Error: Found '//trim(start_st)//' more than once in input file', comm) return end if found_s = .true. end do do loop = 1, settings%num_lines ine = index(settings%in_data(loop), trim(keyword)) if (ine == 0) cycle in = index(settings%in_data(loop), 'end') if (in == 0 .or. in > 1) cycle line_e = loop if (found_e) then call set_error_input(error, & 'w90_readwrite_get_projections: Found '//trim(end_st)//' more than once in input file', comm) return end if found_e = .true. end do if (.not. found_e) then call set_error_input(error, 'w90_readwrite_get_projections: Found '//trim(start_st) & //' but no '//trim(end_st)//' in input file', comm) return end if if (line_e <= line_s) then call set_error_input(error, & 'w90_readwrite_get_projections: '//trim(end_st)//' comes before '//trim(start_st) & //' in input file', comm) return end if block_line_s = line_s block_line_e = line_e dummy = settings%in_data(line_s + 1) if (index(dummy, 'ang') .ne. 0) then if (.not. lcount) settings%in_data(line_s) (1:maxlen) = ' ' line_s = line_s + 1 elseif (index(dummy, 'bohr') .ne. 0) then if (.not. lcount) settings%in_data(line_s) (1:maxlen) = ' ' line_s = line_s + 1 lconvert = .true. elseif (index(dummy, 'random') .ne. 0) then if (.not. lcount) settings%in_data(line_s) (1:maxlen) = ' ' line_s = line_s + 1 if (index(settings%in_data(line_s + 1), end_st) .ne. 0) then lrandom = .true. ! all projections random else lpartrandom = .true. ! only some projections random if (index(settings%in_data(line_s + 1), 'ang') .ne. 0) then if (.not. lcount) settings%in_data(line_s) (1:maxlen) = ' ' line_s = line_s + 1 elseif (index(settings%in_data(line_s + 1), 'bohr') .ne. 0) then if (.not. lcount) settings%in_data(line_s) (1:maxlen) = ' ' line_s = line_s + 1 lconvert = .true. end if end if end if ! skip the begin/end lines line_s = line_s + 1 line_e = line_e - 1 elseif (allocated(settings%entries)) then ! reading from setopt do loop = 1, settings%num_entries if (settings%entries(loop)%keyword == 'projections') then if (settings%entries(loop)%txtdata == 'bohr') lconvert = .true. if (settings%entries(loop)%txtdata == 'random') lrandom = .true. end if end do line_s = 1 line_e = settings%num_entries end if ! reading from input file or entries counter = 0 if (.not. lrandom) then do line = line_s, line_e ang_states = 0 !Assume the default values proj_z_tmp = proj_z_def proj_x_tmp = proj_x_def proj_zona_tmp = proj_zona_def proj_radial_tmp = proj_radial_def if (spinors) then proj_s_qaxis_tmp = proj_s_qaxis_def spn_counter = 2 proj_u_tmp = .true. proj_d_tmp = .true. else spn_counter = 1 end if ! Strip input line of all spaces if (allocated(settings%entries)) then if (settings%entries(line)%keyword /= 'projections') cycle dummy = utility_strip(settings%entries(line)%txtdata) else dummy = utility_strip(settings%in_data(line)) end if if (len(trim(dummy)) == 0) cycle dummy = adjustl(dummy) pos1 = index(dummy, ':') if (pos1 == 0) then call set_error_input(error, & 'w90_wannier90_readwrite_read_projection: malformed projection definition: '//trim(dummy), comm) return end if sites = 0 ctemp = dummy(:pos1 - 1) ! Read the atomic site if (index(ctemp, 'c=') > 0) then sites = -1 ctemp = ctemp(3:) call utility_string_to_coord(ctemp, pos_cart, error, comm) if (allocated(error)) return if (lconvert) pos_cart = pos_cart*bohr call utility_cart_to_frac(pos_cart(:), pos_frac(:), inv_lattice) elseif (index(ctemp, 'f=') > 0) then sites = -1 ctemp = ctemp(3:) call utility_string_to_coord(ctemp, pos_frac, error, comm) if (allocated(error)) return else if (atom_data%num_species == 0) then call set_error_input(error, 'w90_wannier90_readwrite_read_projection: ' & //'Atom centred projection requested but no atoms defined', comm) return end if do loop = 1, atom_data%num_species if (trim(ctemp) == atom_data%label(loop)) then species = loop sites = atom_data%species_num(loop) exit end if if (loop == atom_data%num_species) then call set_error_input(error, 'w90_wannier90_readwrite_read_projection: Atom site not recognised '//trim(ctemp), comm) return end if end do end if dummy = dummy(pos1 + 1:) ! scan for quantisation direction pos1 = index(dummy, '[') if (spinors) then if (pos1 > 0) then ctemp = (dummy(pos1 + 1:)) pos2 = index(ctemp, ']') if (pos2 == 0) then call set_error_input(error, & 'w90_readwrite_get_projections: no closing square bracket for spin quantisation dir', comm) return end if ctemp = ctemp(:pos2 - 1) call utility_string_to_coord(ctemp, proj_s_qaxis_tmp, error, comm) dummy = dummy(:pos1 - 1) ! remove [ ] section end if else if (pos1 > 0) then call set_error_input(error, 'w90_readwrite_get_projections: spin qdir is defined but spinors=.false.', comm) return end if end if ! scan for up or down staring from the end of the string. pos1 = index(dummy, '(', BACK=.true.) ! We need to exclude the case in which we have no spinor specification (u) (d) etc ! But we have an f-orbital specified. if (pos1 > 0) then found_f = .false. ctemp = (dummy(pos1:)) pos2 = index(ctemp, '(x2-y2)') if (pos2 > 0) found_f = .true. pos2 = index(ctemp, '(x2-3y2)') if (pos2 > 0) found_f = .true. pos2 = index(ctemp, '(3x2-y2)') if (pos2 > 0) found_f = .true. if (found_f) pos1 = 0 end if if (spinors) then if (pos1 > 0) then proj_u_tmp = .false.; proj_d_tmp = .false. ctemp = (dummy(pos1 + 1:)) pos2 = index(ctemp, ')') if (pos2 == 0) then call set_error_input(error, 'w90_readwrite_get_projections: no closing bracket for spin', comm) return end if ctemp = ctemp(:pos2 - 1) if (index(ctemp, 'u') > 0) proj_u_tmp = .true. if (index(ctemp, 'd') > 0) proj_d_tmp = .true. if (proj_u_tmp .and. proj_d_tmp) then spn_counter = 2 elseif (.not. proj_u_tmp .and. .not. proj_d_tmp) then call set_error_input(error, 'w90_readwrite_get_projections: found brackets but neither u or d', comm) return else spn_counter = 1 end if dummy = dummy(:pos1 - 1) ! remove ( ) section end if else if (pos1 > 0) then call set_error_input(error, 'w90_readwrite_get_projections: spin is defined but spinors=.false.', comm) return end if end if !Now we know the sites for this line. Get the angular momentum states pos1 = index(dummy, ':') if (pos1 > 0) then ctemp = dummy(:pos1 - 1) else ctemp = dummy end if ctemp2 = ctemp do pos2 = index(ctemp2, ';') if (pos2 == 0) then ctemp3 = ctemp2 else ctemp3 = ctemp2(:pos2 - 1) end if if (index(ctemp3, 'l=') == 1) then mstate = index(ctemp3, ',') if (mstate > 0) then read (ctemp3(3:mstate - 1), *, err=101, end=101) l_tmp else read (ctemp3(3:), *, err=101, end=101) l_tmp end if if (l_tmp < -5 .or. l_tmp > 3) then call set_error_input(error, 'w90_readwrite_get_projections: Incorrect l state requested', comm) return end if if (mstate == 0) then if (l_tmp >= 0) then do loop_m = 1, 2*l_tmp + 1 ang_states(loop_m, l_tmp) = 1 end do elseif (l_tmp == -1) then !sp ang_states(1:2, l_tmp) = 1 elseif (l_tmp == -2) then !sp2 ang_states(1:3, l_tmp) = 1 elseif (l_tmp == -3) then !sp3 ang_states(1:4, l_tmp) = 1 elseif (l_tmp == -4) then !sp3d ang_states(1:5, l_tmp) = 1 elseif (l_tmp == -5) then !sp3d2 ang_states(1:6, l_tmp) = 1 end if else if (index(ctemp3, 'mr=') /= mstate + 1) then call set_error_input(error, 'w90_readwrite_get_projections: Problem reading m state', comm) return end if ctemp4 = ctemp3(mstate + 4:) do pos3 = index(ctemp4, ',') if (pos3 == 0) then ctemp5 = ctemp4 else ctemp5 = ctemp4(:pos3 - 1) end if read (ctemp5(1:), *, err=102, end=102) m_tmp if (l_tmp >= 0) then if ((m_tmp > 2*l_tmp + 1) .or. (m_tmp <= 0)) then call set_error_input(error, 'w90_readwrite_get_projections: m is > l !', comm) return end if elseif (l_tmp == -1 .and. (m_tmp > 2 .or. m_tmp <= 0)) then call set_error_input(error, 'w90_readwrite_get_projections: m has incorrect value (1)', comm) return elseif (l_tmp == -2 .and. (m_tmp > 3 .or. m_tmp <= 0)) then call set_error_input(error, 'w90_readwrite_get_projections: m has incorrect value (2)', comm) return elseif (l_tmp == -3 .and. (m_tmp > 4 .or. m_tmp <= 0)) then call set_error_input(error, 'w90_readwrite_get_projections: m has incorrect value (3)', comm) return elseif (l_tmp == -4 .and. (m_tmp > 5 .or. m_tmp <= 0)) then call set_error_input(error, 'w90_readwrite_get_projections: m has incorrect value (4)', comm) return elseif (l_tmp == -5 .and. (m_tmp > 6 .or. m_tmp <= 0)) then call set_error_input(error, 'w90_readwrite_get_projections: m has incorrect value (5)', comm) return end if ang_states(m_tmp, l_tmp) = 1 if (pos3 == 0) exit ctemp4 = ctemp4(pos3 + 1:) end do end if else do pos3 = index(ctemp3, ',') if (pos3 == 0) then ctemp4 = ctemp3 else ctemp4 = ctemp3(:pos3 - 1) end if read (ctemp4(1:), *, err=106, end=106) m_string select case (trim(adjustl(m_string))) case ('s') ang_states(1, 0) = 1 case ('p') ang_states(1:3, 1) = 1 case ('pz') ang_states(1, 1) = 1 case ('px') ang_states(2, 1) = 1 case ('py') ang_states(3, 1) = 1 case ('d') ang_states(1:5, 2) = 1 case ('dz2') ang_states(1, 2) = 1 case ('dxz') ang_states(2, 2) = 1 case ('dyz') ang_states(3, 2) = 1 case ('dx2-y2') ang_states(4, 2) = 1 case ('dxy') ang_states(5, 2) = 1 case ('f') ang_states(1:7, 3) = 1 case ('fz3') ang_states(1, 3) = 1 case ('fxz2') ang_states(2, 3) = 1 case ('fyz2') ang_states(3, 3) = 1 case ('fxyz') ang_states(4, 3) = 1 case ('fz(x2-y2)') ang_states(5, 3) = 1 case ('fx(x2-3y2)') ang_states(6, 3) = 1 case ('fy(3x2-y2)') ang_states(7, 3) = 1 case ('sp') ang_states(1:2, -1) = 1 case ('sp-1') ang_states(1, -1) = 1 case ('sp-2') ang_states(2, -1) = 1 case ('sp2') ang_states(1:3, -2) = 1 case ('sp2-1') ang_states(1, -2) = 1 case ('sp2-2') ang_states(2, -2) = 1 case ('sp2-3') ang_states(3, -2) = 1 case ('sp3') ang_states(1:4, -3) = 1 case ('sp3-1') ang_states(1, -3) = 1 case ('sp3-2') ang_states(2, -3) = 1 case ('sp3-3') ang_states(3, -3) = 1 case ('sp3-4') ang_states(4, -3) = 1 case ('sp3d') ang_states(1:5, -4) = 1 case ('sp3d-1') ang_states(1, -4) = 1 case ('sp3d-2') ang_states(2, -4) = 1 case ('sp3d-3') ang_states(3, -4) = 1 case ('sp3d-4') ang_states(4, -4) = 1 case ('sp3d-5') ang_states(5, -4) = 1 case ('sp3d2') ang_states(1:6, -5) = 1 case ('sp3d2-1') ang_states(1, -5) = 1 case ('sp3d2-2') ang_states(2, -5) = 1 case ('sp3d2-3') ang_states(3, -5) = 1 case ('sp3d2-4') ang_states(4, -5) = 1 case ('sp3d2-5') ang_states(5, -5) = 1 case ('sp3d2-6') ang_states(6, -5) = 1 case default call set_error_input(error, 'w90_readwrite_get_projections: Problem reading l state '//trim(ctemp3), comm) return end select if (pos3 == 0) exit ctemp3 = ctemp3(pos3 + 1:) end do end if if (pos2 == 0) exit ctemp2 = ctemp2(pos2 + 1:) end do ! check for non-default values if (pos1 > 0) then dummy = dummy(pos1 + 1:) ! z axis pos1 = index(dummy, 'z=') if (pos1 > 0) then ctemp = (dummy(pos1 + 2:)) pos2 = index(ctemp, ':') if (pos2 > 0) ctemp = ctemp(:pos2 - 1) call utility_string_to_coord(ctemp, proj_z_tmp, error, comm) if (allocated(error)) return end if ! x axis pos1 = index(dummy, 'x=') if (pos1 > 0) then ctemp = (dummy(pos1 + 2:)) pos2 = index(ctemp, ':') if (pos2 > 0) ctemp = ctemp(:pos2 - 1) call utility_string_to_coord(ctemp, proj_x_tmp, error, comm) if (allocated(error)) return end if ! diffusivity of orbital pos1 = index(dummy, 'zona=') if (pos1 > 0) then ctemp = (dummy(pos1 + 5:)) pos2 = index(ctemp, ':') if (pos2 > 0) ctemp = ctemp(:pos2 - 1) read (ctemp, *, err=104, end=104) proj_zona_tmp end if ! nodes for the radial part pos1 = index(dummy, 'r=') if (pos1 > 0) then ctemp = (dummy(pos1 + 2:)) pos2 = index(ctemp, ':') if (pos2 > 0) ctemp = ctemp(:pos2 - 1) read (ctemp, *, err=105, end=105) proj_radial_tmp end if end if ! if (sites == -1) then ! if (counter + spn_counter*sum(ang_states) > num_proj) & ! call io_error('w90_readwrite_get_projections: too many projections defined') ! else ! if (counter + spn_counter*sites*sum(ang_states) > num_proj) & ! call io_error('w90_readwrite_get_projections: too many projections defined') ! end if if (sites == -1) then do loop_l = min_l, max_l do loop_m = min_m, max_m if (ang_states(loop_m, loop_l) == 1) then do loop_s = 1, spn_counter counter = counter + 1 if (lcount) cycle input_proj(counter)%site(:) = pos_frac input_proj(counter)%l = loop_l input_proj(counter)%m = loop_m input_proj(counter)%z(:) = proj_z_tmp input_proj(counter)%x(:) = proj_x_tmp input_proj(counter)%radial = proj_radial_tmp input_proj(counter)%zona = proj_zona_tmp if (spinors) then if (spn_counter == 1) then if (proj_u_tmp) input_proj(counter)%s = 1 if (proj_d_tmp) input_proj(counter)%s = -1 else if (loop_s == 1) input_proj(counter)%s = 1 if (loop_s == 2) input_proj(counter)%s = -1 end if input_proj(counter)%s_qaxis(:) = proj_s_qaxis_tmp end if end do end if end do end do else do loop_sites = 1, sites do loop_l = min_l, max_l do loop_m = min_m, max_m if (ang_states(loop_m, loop_l) == 1) then do loop_s = 1, spn_counter counter = counter + 1 if (lcount) cycle call utility_cart_to_frac(atom_data%pos_cart(:, loop_sites, species), & pos_frac, inv_lattice) input_proj(counter)%site(:) = pos_frac(:) input_proj(counter)%l = loop_l input_proj(counter)%m = loop_m input_proj(counter)%z(:) = proj_z_tmp input_proj(counter)%x(:) = proj_x_tmp input_proj(counter)%radial = proj_radial_tmp input_proj(counter)%zona = proj_zona_tmp if (spinors) then if (spn_counter == 1) then if (proj_u_tmp) input_proj(counter)%s = 1 if (proj_d_tmp) input_proj(counter)%s = -1 else if (loop_s == 1) input_proj(counter)%s = 1 if (loop_s == 2) input_proj(counter)%s = -1 end if input_proj(counter)%s_qaxis(:) = proj_s_qaxis_tmp end if end do end if end do end do end do end if end do !end loop over projection block ! check there are enough projections and add random projections if required if (.not. lpartrandom) then if (counter .lt. num_wann) then call set_error_input(error, 'w90_readwrite_get_projections: too few projection functions defined', comm) return end if end if end if ! .not. lrandom if (lcount) then if (counter .lt. num_wann) then num_proj = num_wann else num_proj = counter end if return end if if (lpartrandom .or. lrandom) then call random_seed() ! comment out this line for reproducible random positions! num_proj = num_wann ! input_proj is allocated 1:num_wann when 'random' do loop = counter + 1, num_wann call random_number(input_proj(loop)%site(:)) input_proj(loop)%l = 0 input_proj(loop)%m = 1 input_proj(loop)%z(:) = proj_z_def input_proj(loop)%x(:) = proj_x_def input_proj(loop)%zona = proj_zona_def input_proj(loop)%radial = proj_radial_def if (spinors) then if (modulo(loop, 2) == 1) then input_proj(loop)%s = 1 else input_proj(loop)%s = -1 end if input_proj(loop)%s_qaxis(1) = 0.0_dp input_proj(loop)%s_qaxis(2) = 0.0_dp input_proj(loop)%s_qaxis(3) = 1.0_dp end if end do end if ! I shouldn't get here, but just in case if (.not. lcount .and. allocated(settings%in_data)) then settings%in_data(block_line_s:block_line_e) (1:maxlen) = ' ' end if !~ ! Check !~ do loop=1,num_wann !~ if ( abs(sum(proj_z(:,loop)*proj_x(:,loop))).gt.1.0e-6_dp ) then !~ write(stdout,*) ' Projection:',loop !~ call io_error(' Error in projections: z and x axes are not orthogonal') !~ endif !~ enddo ! Normalise z-axis and x-axis and check/fix orthogonality do loop = 1, num_proj znorm = sqrt(sum(input_proj(loop)%z(:)*input_proj(loop)%z(:))) xnorm = sqrt(sum(input_proj(loop)%x(:)*input_proj(loop)%x(:))) input_proj(loop)%z(:) = input_proj(loop)%z(:)/znorm ! normalise z input_proj(loop)%x(:) = input_proj(loop)%x(:)/xnorm ! normalise x cosphi = sum(input_proj(loop)%z(:)*input_proj(loop)%x(:)) ! Check whether z-axis and z-axis are orthogonal if (abs(cosphi) .gt. eps6) then ! Special case of circularly symmetric projections (pz, dz2, fz3) ! just choose an x-axis that is perpendicular to the given z-axis if ((input_proj(loop)%l .ge. 0) .and. (input_proj(loop)%m .eq. 1)) then proj_x_tmp(:) = input_proj(loop)%x(:) ! copy of original x-axis call random_seed() call random_number(proj_z_tmp(:)) ! random vector ! calculate new x-axis as the cross (vector) product of random vector with z-axis input_proj(loop)%x(1) = proj_z_tmp(2)*input_proj(loop)%z(3) & - proj_z_tmp(3)*input_proj(loop)%z(2) input_proj(loop)%x(2) = proj_z_tmp(3)*input_proj(loop)%z(1) & - proj_z_tmp(1)*input_proj(loop)%z(3) input_proj(loop)%x(3) = proj_z_tmp(1)*input_proj(loop)%z(2) & - proj_z_tmp(2)*input_proj(loop)%z(1) xnorm_new = sqrt(sum(input_proj(loop)%x(:)*input_proj(loop)%x(:))) input_proj(loop)%x(:) = input_proj(loop)%x(:)/xnorm_new ! normalise goto 555 end if ! If projection axes non-orthogonal enough, then ! user may have made a mistake and should check if (abs(cosphi) .gt. eps2) then write (stdout, *) ' Projection:', loop call set_error_input(error, ' Error in projections: z and x axes are not orthogonal', comm) return end if ! If projection axes are "reasonably orthogonal", project x-axis ! onto plane perpendicular to z-axis to make them more so sinphi = sqrt(1 - cosphi*cosphi) proj_x_tmp(:) = input_proj(loop)%x(:) ! copy of original x-axis ! calculate new x-axis: ! x = z \cross (x_tmp \cross z) / sinphi = ( x_tmp - z(z.x_tmp) ) / sinphi input_proj(loop)%x(:) = (proj_x_tmp(:) - cosphi*input_proj(loop)%z(:))/sinphi ! Final check 555 cosphi_new = sum(input_proj(loop)%z(:)*input_proj(loop)%x(:)) if (abs(cosphi_new) .gt. eps6) then write (stdout, *) ' Projection:' call set_error_input(error, ' Error: z and x axes are still not orthogonal after projection', comm) return end if end if end do return 101 call set_error_input(error, 'w90_readwrite_get_projections: Problem reading l state into integer '//trim(ctemp3), comm) return 102 call set_error_input(error, 'w90_readwrite_get_projections: Problem reading m state into integer '//trim(ctemp3), comm) return 104 call set_error_input(error, 'w90_readwrite_get_projections: Problem reading zona into real '//trim(ctemp), comm) return 105 call set_error_input(error, 'w90_readwrite_get_projections: Problem reading radial state into integer '//trim(ctemp), comm) return 106 call set_error_input(error, 'w90_readwrite_get_projections: Problem reading m state into string '//trim(ctemp3), comm) return end subroutine w90_readwrite_get_projections