1program distances 2! 3!.... Convert a typical PDB XYZ file to a list of internuclear distances 4! 5!.... Program runs as a filter (old-school UNIX), reading from stdin 6!.... and writing to stdout, and is invoked as 7!.... distances [--min lower_bound] [--max upper_bound] < XYZinputfile > Distancesfile 8!.... where the optional arguments screen out internuclear 9!.... distances greater than upper_bound and less than lower_bound. 10! 11!.... If --min is omitted, lower_bound is taken as zero. 12!.... If --max is omitted, upper_bound is taken as (effectively) infinity. 13! 14!.... The less explicit, old-school form 15!.... distances [-b lower_bound] [-t upper_bound] < XYZinputfile > Distancesfile 16!.... will also work. 17! 18!.... Assumptions: 19!.... All atoms of a given atomic number may not be grouped together. 20!.... There are at most 99 atoms of a given type, except for 21!.... those with a one-letter atomic symbol. 22!.... No more than 999 atoms total. 23!.... Highest atomic number handled is 112 (Cn) 24! 25 26 use io_channels 27 use periodic_table 28 29 implicit none 30 31 integer, parameter :: MAX_LENGTH=80 32 character(LEN=MAX_LENGTH) :: input_line 33 character(LEN=MAX_LENGTH) :: output_line(999,MAX_ELTS) 34 character(LEN=MAX_LENGTH) :: temp_line(999,MAX_ELTS) 35 integer :: atoms_of_this_type(MAX_ELTS) 36 37 integer :: num_atoms, num_types 38 39 integer :: arg_count 40 41 character(LEN=2) :: symbol 42 integer :: atomic_number 43 integer :: total_atoms 44 real (kind(0.d0)) :: x, y, z 45 real (kind(0.d0)) :: screening_threshold, baseline_threshold 46 real (kind(0.d0)) :: x_vec(999), y_vec(999), z_vec(999) 47 real (kind(0.d0)) :: dist 48 49 character(LEN=4) :: name_vec(999) 50 51 logical :: com_transform 52 53 54 55 real (kind(0.d0)), parameter :: angstrom_to_a0=0.5291772083D0 56 57 character(LEN=2) :: current_symbol 58 integer :: current_count 59 logical :: one_symbol 60 character(LEN=5) :: charge(MAX_ELTS) 61 character(LEN=5) :: count(MAX_ELTS) 62 63 integer :: i, j, k 64 65! Check for command-line argument for screening threshold 66 67 screening_threshold = 1.d8 68 baseline_threshold = -0.1d0 69 if (command_argument_count() .eq. 2) then 70 call get_command_argument(1,input_line) 71 if (input_line(1:2) .eq. '-t' & 72 .or. input_line(1:5) .eq. '--max') then 73 call get_command_argument(2,input_line) 74 read(input_line,*) screening_threshold 75 elseif (input_line(1:2) .eq. '-b' & 76 .or. input_line(1:5) .eq. '--min') then 77 call get_command_argument(2,input_line) 78 read(input_line,*) baseline_threshold 79 endif 80 endif 81 82 if (command_argument_count() .eq. 4) then 83 call get_command_argument(1,input_line) 84 if (input_line(1:2) .eq. '-t' & 85 .or. input_line(1:5) .eq. '--max') then 86 call get_command_argument(2,input_line) 87 read(input_line,*) screening_threshold 88 elseif (input_line(1:2) .eq. '-b' & 89 .or. input_line(1:5) .eq. '--min') then 90 call get_command_argument(2,input_line) 91 read(input_line,*) baseline_threshold 92 endif 93 94 call get_command_argument(3,input_line) 95 if (input_line(1:2) .eq. '-t' & 96 .or. input_line(1:5) .eq. '--max') then 97 call get_command_argument(4,input_line) 98 read(input_line,*) screening_threshold 99 elseif (input_line(1:2) .eq. '-b' & 100 .or. input_line(1:5) .eq. '--min') then 101 call get_command_argument(4,input_line) 102 read(input_line,*) baseline_threshold 103 endif 104 105 endif 106 107 atoms_of_this_type = 0 108 109 read(luinp,*) num_atoms 110 read(luinp,'(a)') input_line 111 112 num_types = 0 113 114 current_symbol = ' ' 115 current_count = 0 116 117 do i = 1,num_atoms 118 read(luinp,'(a)') input_line 119! write(luerr,'(a)') input_line 120 121 symbol = input_line(1:2) 122 123 if (symbol .eq. current_symbol) then 124 current_count = current_count + 1 125 atoms_of_this_type(atomic_number) = current_count 126 127 read(input_line(3:MAX_LENGTH),*) x, y, z 128 129 x = x/angstrom_to_a0 130 y = y/angstrom_to_a0 131 z = z/angstrom_to_a0 132 133 if (one_symbol) then 134 write(output_line(current_count,atomic_number),'(a1,i3.3,3(4x,f14.6))') & 135 trim(current_symbol), current_count, x, y, z 136 else 137 write(output_line(current_count,atomic_number),'(a2,i2.2,3(4x,f14.6))') & 138 current_symbol, current_count, x, y, z 139 endif 140 141 else 142 143!.... Different element: write out any preceding block before locating it. 144 145 if (current_count .ne. 0) then 146 select case(atomic_number) 147 case(1:9) 148 write(charge(atomic_number),'(f3.1,a)') real(atomic_number)+0.001, ' ' 149 case(10:99) 150 write(charge(atomic_number),'(f4.1,a)') real(atomic_number)+0.001, ' ' 151 case(100:) 152 write(charge(atomic_number),'(f5.1,a)') real(atomic_number)+0.001 153 end select 154 155 select case(current_count) 156 case(1:9) 157 write(count(atomic_number),'(i1)') current_count 158 case(10:99) 159 write(count(atomic_number),'(i2)') current_count 160 case(100:) 161 write(count(atomic_number),'(i3)') current_count 162 end select 163 164 endif 165 166!.... Locate new element 167 do j = 1,MAX_ELTS 168 if (symbol .ne. elements(j)) cycle 169 170 current_symbol = symbol 171 one_symbol = len(trim(symbol)) .eq. 1 172 atomic_number = j 173 if (atoms_of_this_type(atomic_number) .eq. 0) num_types = num_types + 1 174 current_count = atoms_of_this_type(atomic_number) + 1 175 atoms_of_this_type(atomic_number) = current_count 176 177 read(input_line(3:MAX_LENGTH),*) x, y, z 178 179 x = x/angstrom_to_a0 180 y = y/angstrom_to_a0 181 z = z/angstrom_to_a0 182 183 if (one_symbol) then 184 write(output_line(current_count,atomic_number),'(a1,i3.3,3(4x,f14.6))') & 185 trim(current_symbol), current_count, x, y, z 186 else 187 write(output_line(current_count,atomic_number),'(a2,i2.2,3(4x,f14.6))') & 188 current_symbol, current_count, x, y, z 189 endif 190 191 exit 192 enddo 193 194 endif 195 196 enddo 197 198!.... Write last block 199 200 if (current_count .ne. 0) then 201 select case(atomic_number) 202 case(1:9) 203 write(charge(atomic_number),'(f3.1,a)') real(atomic_number)+0.001, ' ' 204 case(10:99) 205 write(charge(atomic_number),'(f4.1,a)') real(atomic_number)+0.001, ' ' 206 case(100:) 207 write(charge(atomic_number),'(f5.1,a)') real(atomic_number)+0.001 208 end select 209 210 select case(current_count) 211 case(1:9) 212 write(count(atomic_number),'(i1)') current_count 213 case(10:99) 214 write(count(atomic_number),'(i2)') current_count 215 case(100:) 216 write(count(atomic_number),'(i3)') current_count 217 end select 218 219 endif 220 221! Now generate internuclear distances and write out 222 223 total_atoms = 0 224 225 do atomic_number = MAX_ELTS,1,-1 226 227 if (atoms_of_this_type(atomic_number) .gt. 0) then 228 229 230 do j = 1,atoms_of_this_type(atomic_number) 231 232 total_atoms = total_atoms + 1 233 234 read(output_line(j,atomic_number),'(a4,3(4x,f14.6))') & 235 name_vec(total_atoms), x_vec(total_atoms), & 236 y_vec(total_atoms), z_vec(total_atoms) 237 238 enddo 239 240 endif 241 242 enddo 243 244 write(lupri,'(/,a,/)') 'Table of internuclear distances in Angstrom' 245 write(lupri,'(4x,a4,4x,a4,14x,a8,/)') 'Atom', 'Atom', 'Distance' 246 247 do j = 2,total_atoms 248 249 do k = 1,j-1 250 251 dist = angstrom_to_a0 & 252 *sqrt((x_vec(k) - x_vec(j))**2 & 253 + (y_vec(k) - y_vec(j))**2 & 254 + (z_vec(k) - z_vec(j))**2) 255 256 if (dist .lt. screening_threshold .and. dist .gt. baseline_threshold) & 257 write(lupri,'(4x,a4,4x,a4,8x,f14.6)') & 258 name_vec(j), name_vec(k), dist 259 260 enddo 261 262 enddo 263 264end program distances 265 266 267 268