First off, your problem with  $rest1[$i] might be due to this:

$i=0;
while ($line=<SOLVFILE>)
{
    chop($line);  # use "chomp" here, not "chop"
    ...
}
[download]

BTW, what value do you see when you print "Number of atoms in solvent file is $Nsolv \n" ?

As for general improvements, many people prefer to feed file names to the script by putting them on the command line with the name of the script:

  my_script.pl  file1  file2
[download]

Personally, I would prefer to structure the input data to cut down on the duplicated code -- your two input files seem to be organized the same way, so you could read both files into a single set of hashes, keyed by file name; each hash element would be an array of values to be used from the file (atom names, coords, radii). This means you can loop over the file names, and the code that reads and structures the data is needed in only on place in the script.

(Oh, as I worked out the concept for myself, it looked to me like you have an error in the "sqrt" call -- you use the "z" values twice, and don't use the "y" values at all. Restructuring the code and data storage to reduce the amount of typing will help to avoid such mistakes.)

Finally, proper indentation really does help. Here's one approach -- it parses properly, but I don't have data to test it on:

#!/usr/bin/perl -w

use strict; # painful at first, maybe, but worth it

# using command line args is good (and easy):

die "Usage: $0 solvent.file structure.file\n" unless
    ( @ARGV == 2 and -f $ARGV[0] and -f $ARGV[1] );

open( PDB, ">Solvent.pdb" ) or die "...:$!\n";
open( TXT, ">OverL.txt" ) or die "...:$!\n";

my %input;
my %atmnb;   # declare storage near where you need it
my %coord;
my %natms;

for my $file ( @ARGV ) {
    open( IN, $file ) or die "open failed on $file: $!\n";

    while (<IN>) { # use a slice on split to get just the needed data:
        my ($atb,$x,$y,$z,$rad) = (split)[1,5,6,7,9];
        push @{$input{$file}}, $_;  # keep full input line
        push @{$atmnb{$file}}, $atb;
        push @{$coord{$file}}, [ $x, $y, $z, $rad ];
    }
    close IN;
    print "Number of atoms in $file is $.\n";
    $natms{$file} = $.;
}

my ($f1,$f2) = @ARGV;

for my $i ( 0 .. $natms{$f1}-1 ) {
    for my $j ( 0 .. $natms{$f2}-1 ) {
# use map to compute the x,y,z diffs just once each:
        my @difs = map { $coord{$f1}[$i][$_] - $coord{$f2}[$j][$_] } (
+ 0, 1, 2 );
        my $dist = sqrt( $difs[0] * $difs[0] +
                         $difs[1] * $difs[1] +
                         $difs[2] * $difs[2] );
        my $rsum = $coord{$f1}[$i][3] + $coord{$f2}[$j][3];
        if ( $dist < $rsum ) {
            print TXT "Overlap for $atmnb{$f1}[$i] and $atmnb{$f2}[$j]
+\n";
        }
        else {
            print PDB "$input{$f1}[$i]";
        }
    }
}
[download]

Another update: regarding optimization, I don't know how much it would save, but in that N1 * N2 nested loop for measuring the distances, you could extract the "f1" data into simpler storage, and save the work of having to dereference the "f1" data structure on every iteration of the inner loop -- ie:

for my $i ( 0 .. $natms{$f1}-1 ) {

    my $xyz = @{$coord{$f1}[$i]};
    my $rad = pop @xyz;

    for my $j ( 0 .. $natms{$f2}-1 ) {
# use map to compute the x,y,z diffs just once each:
        my @difs = map { $xyz[$_] - $coord{$f2}[$j][$_] } ( 0, 1, 2 );
        my $dist = sqrt( $difs[0] * $difs[0] +
                         $difs[1] * $difs[1] +
                         $difs[2] * $difs[2] );
        my $rsum = $rad + $coord{$f2}[$j][3];
        if ( $dist < $rsum ) {
            print TXT "Overlap for $atmnb{$f1}[$i] and $atmnb{$f2}[$j]
+\n";
        }
        else {
            print PDB $input{$f1}[$i];
        }
    }
}
[download]

If you're on a system that uses "\r\n" for line termination (i.e. windows), "chop" will only take away the "\n", leaving "\r", whereas "chomp" will take away both.

Wrong. If you are Win32, then "\r\n" is used for line termination in files and the record you read in (unless you did binmode) will end with just "\n". $/ starts out as "\n" everywhere. See Re^4: Line Feeds (rumor control) for more on this.

The reason chomp is preferred to chop for this is because chop will remove the last character of the line even if it isn't a newline.

                - tye

somewhat optimised - I like to think of it as a good mix between obvious and efficient.

 #!/usr/bin/perl -w

use strict;

my @solvents;
my @structs;

my @fields = qw(atom number name resname resnumber xpos ypos zpos q ra
+dius);

slurp_file(\@solvents, "solv.txt");
slurp_file(\@structs,  "struct.txt");

print "Atoms in solvent file is ",   scalar(@solvents), "\n";
print "Atoms in structure file is ", scalar(@structs),  "\n";

my @desBooli = map {0} @solvents; # could be done other ways

open(test2, ">OverL.txt")
  or die ("Couldn't open file OverL.txt :: $!\n");

for (my $i = 0; $i < @solvents; $i++) {
  for (my $j = 0; $j < @structs; $j++) {
    # I am sure there is a better way
    my $dist = sqrt(($solvents[$i]->{xpos} - $structs[$j]->{xpos})**2+
                    ($solvents[$i]->{ypos} - $structs[$j]->{ypos})**2+
                    ($solvents[$i]->{zpos} - $structs[$j]->{zpos})**2)
+;
    my $radSum = $solvents[$i]->{radius} + $structs[$j]->{radius};
    if($dist < $radSum) {
      print test2 "$i - Overlap for atoms solvent $solvents[$i]->{resn
+umber} and structure $structs[$j]->{resnumber}\n";
      $desBooli[$i] = 1;
    }
  }
}
close test2;

open(test1, ">Solvent.pdb")
  or die ("Couldn't open file Solvent.pdb :: $!\n");
for (my $i = 0; $i < @solvents; $i++) {
    if ($desBooli[$i] == 0){ # no overlap
       print test1 dump_solvent($solvents[$i]), "\n";
    }
}
close test1;

exit;

sub slurp_file {
  my ($array, $name) = @_;

  open(FILE, $name) or die ("Couldn't open file $name:: $!\n");
  while (<FILE>) {
    chomp; # trim line endings
    next if /^\s*$/; # strip blank lines
    my %hash = ();
    @hash{@fields} = split /\s+/; # hash slice
    push @{$array}, \%hash;
  }
  use Data::Dumper;
  close FILE;
}

sub dump_solvent {
  my ($solvent) = @_;

  join ' ', @{$solvent}{@fields};
}
[download]

I'm a bit reluctant to go into this because implementing this would be a rather complex programming challenge. But I'll explain it as best I can in hopes that someone will find it helpful, educational, or at least interesting. (:

When I previously had to find overlaps between two sets of shapes, I certainly couldn't wait to compare against every single shape in the database. So I had to come up with an index on the table that would allow me to quickly eliminate almost all of the items in the database as not being nearly close enough to have to worry about.

The problem with a database index is that it is one-dimensional.

I had previously had to deal with finding the nearest point (not overlapping shapes) and for that I really needed a two-dimensional key. But after some work I figured out that I could do pretty good by sorting the points into parallel "bands" that would define a linear (one-dimensional) ordering of all the points and allow me to build an index.

So, if I made the bands W units thick and horizontal and wanted to find all point within 2*W units of some point, then I might end up doing a keyed search for records falling into the following 5 sections of bands:

----------------------------
             [ ]
----------------------------
        [           ]
----------------------------
    [         .          ]
----------------------------
       [             ]
----------------------------
            [   ]
----------------------------
[download]

For finding overlapping shapes, things get a bit more complicated. You can key each shape based on it's center. But you'd have to be able to assert that each shape is no larger than some maximum size or else you couldn't restrict your search.

Based on your problem description, I think you might be able to just figure out your maximum radius and index your points into "shafts" and have a pretty fast solution. That is, for each point in the smaller set (or in the larger set if this set is used over and over with few changes against a bunch of different smaller sets) you build a key from the following values (in the order given):

    int(X/W), int(Y/W), Z
[download]

Then, if you need to check a sphere centered at (X1,Y1,Z1) and with radius R1, you need to look for all spheres in the database where their center is within R1+M of (X1,Y1,Z1). So loop over all int(X/W) between X1-R1-M and X1+R1+M:

for my $xBand (  int(($X1-$R1-$M)/$W) .. int(($X1-$R1-$M)/$W)  ) {
[download]

  for my $yBand (  int(($Y1-$R1-$M)/W) .. int(($Y1-$R1-$M)/$W)  ) {
[download]

Now you are looping over all shafts that contain centers within R1+M of the current center of interest. So, for each shaft, look for centers in the range. We'll again not optimize this for the sake of simplicity:

    my $start= BuildKey( $xBand, $yBand, $Z1-$R1-$M );
    my $end=   BuildKey( $xBand, $yBand, $Z1+$R1+$M );
    ... "select * from shapes where key between $start and $end"
    my $record;
    while(  $record= Fetch( ... )   ) {
        if(  Overlap( $record, ... )  ) {
            Report( ... );
            return 1;
        }
    }
  }
}
return 0;
[download]

But, what if you can't pick a reasonable maximum radius? Well, you just segregate shapes based on what range their radius is in and add a loop around all of the above code that iterates over each radius range.

For example, you might realize that most of your radii are <= 2, a fraction of them are between 2 and 10, and a very few are over 10. In this case you might build a key as follows:

my @bandWidth= ( 3, 15, 0 );

sub BuildKey {
    my( $radius, $x, $y, $z )= @_;
    my $sizeRange=
        $radius <= 2 ? 0
      : $radius <= 10 ? 1
      : 2;
    my $width= $bandWidth[$sizeRange];
    if(  0 == $width  ) {
        return ConcatForKey( $sizeRange, 0, 0, 0 );
    }
    return ConcatForKey(
        $sizeRange, int($x/$width), int($y/$width), $z );
}
[download]

    real        x
    real        y
    real        z
    real        radius
    int         size
    int         xBand
    int         yBand
    index       size,xBand,yBand,z
[download]

Or you chould use a DBM file in which case you (pretty much) need to express each key as a string (of bytes) such that normal "ASCIIbetical" ordering (as done by cmp) sorts properly.

Assuming that your X and Y coordinates never go beyond -2147483647*$Width nor 2147483647*$Width, then you can put big-endian longs into your key string and get a proper sort order. However, you have to shift the values so that negative numbers sort before positive numbers. We'll even put the Z values into "bands" so we can use the same strategy for putting them in the key that we do for the X and Y band numbers:

sub BuildKey {
    my( $radius, $x, $y, $z )= @_;
    my $sizeRange=
        $radius <= 2 ? 0
      : $radius <= 10 ? 1
      : 2;
    my $width= $bandWidth[$sizeRange];
    if(  0 == $width  ) {
        return pack "c", $sizeRange;
    }
    for my $c (  $x, $y, $z  ) {
        $c= int($c/$width);
        $c ^= 0x8000_0000;
    }
    return return pack "cNNN", $sizeRange, $x, $y, $z;
}
[download]

Note that this assumes 4-byte integers. If you have a version of Perl built with 8-byte integer support, then you'll probably have to change the ^= line. If you want to use 8-byte integers in your keys, then you'll also have to adjust the pack expression.

I don't see anything specifically wrong with line 58 and the final for loop. However I did find a serious typo, in your code to find the distance between the two atoms, you accidentally used the Z co-ordinate twice instead of Y. I'm not able to help you optimize the algorithm, but I did clean up the code so it is much more readable (and runs under strict). Allthough I did use 'constant' so all of the caveats for it apply. Be careful of PM wrapping this code. Oh and since you only provided one example line I tested this against some randomly generated data which may not have been correct :).

#!/usr/bin/perl -w
use strict;

use constant {
            RAD => 9,
            X => 5,
            Y => 6,
            Z => 7,
            NB => 1,
            REST => 10,
};

my $file;
print "Please type the filename with the solvent: ";
chomp($file = <STDIN>);
die "File \"$file\" doesn\'t seem to exist!!\n" unless -e $file;
open(SOLVFILE, $file) or die("Couldn't open file  $file\n");

print "Please type the filename with the structure: \n";
chomp($file = <STDIN>);
die "File \"$file\" doesn\'t seem to exist!!\n" unless -e $file;
open(STRUCTFILE, $file) or die("Couldn't open file  $file\n");

open(PDB, ">Solvent.pdb") or die "$!";
open(OVERL, ">OverL.txt") or die "$!";

my @solv;
my $line;
while ($line=<SOLVFILE>)
{
    chomp $line;
    push @solv, [split /\s+/, $line];
    push @{$solv[-1]}, $line;
}
close SOLVFILE;
print "Number of atoms in solvent file is ",$#solv+1," \n";

my @struct;
while ($line=<STRUCTFILE>)
{
    chomp($line);
    push @struct, [split /\s+/, $line];
    push @{$struct[-1]}, $line;
}
close STRUCTFILE;
print "Number of atoms in structure file is ",$#struct+1," \n";

for (my $i = 0; $i <= $#solv; $i++)
{
    my $v = $solv[$i];
    for my $t (@struct)
    {
        
        if (    # Sum of radius is greater than distance
                # between the atoms
            ($v->[RAD] + $t->[RAD])
            > 
            sqrt(
              ($v->[X] - $t->[X])**2
            + ($v->[Y] - $t->[Y])**2
            + ($v->[Z] - $t->[Z])**2
            ))
        {
            print OVERL "Overlap for atoms ",
                   $v->[NB]," and ",$t->[NB],"\n";
        } else {
                print PDB $v->[REST],"\n";
        }
    }
}
close OVERL;
close PDB;
exit;
[download]

HTH

tedrek

Update: I realized the whole @overlap functionality could be moved into the main loop, so I did :). Also you could change the outer for loop to to a while(<SOLV>) and not need to read SOLV into memory.

Update2: After seeing graff's code I realized the print to PDB was supposed to be in a else block.