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