Расстояние между одной точкой и всеми остальными в файле PDB

Не уверен, что я четко понимаю, что вы хотите, но как насчет:

#!/usr/local/bin/perl 
use strict;
use warnings;

my (@refer, @points);
my $part = 0;
while (my $line = <DATA>) { 
    chomp($line);
    if ($line =~ /^TER/) {
        $part++;
        next;
    }
    my @array = (split (/\s+/, $line))[5, 6, 7]; 
    if ($part == 0) {
        push @refer, [ @array ]; 
    } else {
        push @points, [ @array ]; 
    }
} 
my %max = (val=>0, x=>0, y=>0, z=>0);
foreach my $ref(@refer) {
    my ($x1, $y1, $z1) = @{$ref};
    foreach my $atom(@points) {
        my ($x, $y, $z) = @{$atom};
        my $dist = sqrt( ($x-$x1)**2 + ($y-$y1)**2 + ($z-$z1)**2 );
        if ($dist > $max{val}) {
            $max{val} = $dist;
            $max{x} = $x;
            $max{y} = $y;
            $max{z} = $z;
        }
    }
}
print "max is $max{val}; coord: x=$max{x}, y=$max{y}, z=$max{z}\n";

__DATA__
ATOM   1279 C    ALA    81      -1.925 -11.270   1.404
ATOM   1280 O    ALA    81      -0.279   9.355  15.557
ATOM   1281 OXT  ALA    81      -2.188  10.341  15.346
TER   
ATOM   1282 N    THR    82      29.632   5.205   5.525
ATOM   1283 H1   THR    82      30.175   4.389   5.768
ATOM   1284 H2   THR    82      28.816   4.910   5.008

вывод:

max is 35.9813670807545; coord: x=30.175, y=4.389, z=5.768

Основная проблема здесь - чтение данных. Во-первых, обратите внимание, что нельзя использовать разделение с текстовыми файлами PDB, поскольку поля определяются по положению, а не по разделителям. См. Описание файла координат (формат PDB).

Чтобы разделить запись ATOM для различных цепочек полимеров, вы можете начать с упрощенной версии, например

my $iblock = 0;
my @atoms = ();
while (my $line = <IN>) { 
   chomp($line);

   # Switch blocks at TER lines
   if ($line =~ /^TER/) {
      $iblock++;

   # Read ATOM lines
   } elsif ($line =~ m/^ATOM/) {
      my @xyz = (substr($line,7-1,9),substr($line,16-1,9),substr($line,25-1,9)); 
      printf "Block %d: atom at (%s)\n",$iblock,join (",",@xyz); 
      push @{$atoms[$iblock]},\@xyz; 

   # Parse additional line types (if needed)
   } else {
      ...
   }
} 

Затем следует цикл по всем парам координат из разных блоков, структурированный следующим образом:

# 1st block
for my $iblock1 (0..$#atoms) {

   # 2nd block
   for my $iblock2 ($iblock1+1..$#atoms) {

      # Compare all pairs of atoms
      ... 
      my $xyz1 (@{$atoms[$iblock1]}) {
         for my $xyz2 (@{$atoms[$iblock2]}) {
            # Calculate distance and compare with $max_dist
            ...
         }
      }
      # Print the maximal distance between these two blocks
      ...
   }
}

Конечно, код может быть более общим, если используется более сложная структура данных или один из доступных парсеров PDB, например Bioperl.

При правильной инкапсуляции это довольно просто и требует незначительных изменений вашего кода.

ETA: добавлено решение с фиксированной шириной, которое у меня было под рукой. Вероятно, было бы лучше прочитать все поля вместо того, чтобы отбрасывать первые 31 символ, а затем вернуть их все в хеш-ссылке. Таким образом, вы можете обрабатывать все строки с помощью одной и той же подпрограммы и просто переключаться между частями, когда первое поле оказывается TER. Вам должно быть легко экстраполировать это из данного кода.

Вы заметите, что ссылочные значения считываются с помощью цикла, потому что нам нужно прервать цикл в точке останова. Остальные значения записываются с помощью оператора map. Затем мы просто передаем данные подпрограмме, которую создали из вашего исходного кода (с некоторыми улучшениями). Я использовал те же имена для лексических переменных, чтобы облегчить чтение кода.

use strict;
use warnings;

my @points;
while (<DATA>) {
    last if /^TER$/;
    push @points, getpoints($_);
}
my @ref = map getpoints($_), <DATA>;

for my $p (@points) {
    getcoords($p, \@ref);
}

sub getpoints {
    my $line = shift;
    my @data = unpack "A31 A8 A8 A8", $line;
    shift @data;
    return \@data;
}
sub getcoords {
    my ($p, $ref) = @_;
    my ($p1,$p2,$p3) = @$p;
    my $max=0;
    my ($x,$y,$z);
    for my $aref ( @$ref ) {
        my ( $x1, $y1, $z1 ) = @$aref;  
        my $dist = sqrt(
            ($x1-$p1)**2 +
            ($y1-$p2)**2 +
            ($z1-$p3)**2
        ); 
        print "distance from ($p1 $p2 $p3) to ( $x1, $y1, $z1 ) is $dist\n"; 

        if ( $dist > $max ) {
            $max = $dist;
            $x=$x1;
            $y=$y1;
            $z=$z1; 
        }
    }
    print "maximum value is : $max\n";
    print "co ordinates are : $x $y $z\n";
}

__DATA__
ATOM   1279 C    ALA    81      -1.925 -11.270   1.404
ATOM   1280 O    ALA    81      -0.279   9.355  15.557
ATOM   1281 OXT  ALA    81      -2.188  10.341  15.346
TER
ATOM   1282 N    THR    82      29.632   5.205   5.525
ATOM   1283 H1   THR    82      30.175   4.389   5.768
ATOM   1284 H2   THR    82      28.816   4.910   5.008