#!/usr/bin/perl
#
#  program:     lammps2pdb.pl
#  author:      Pieter J. in 't Veld
#               pjintve@sandia.gov, veld@verizon.net
#  date:        September 3-4, 2004, April 23-25, 2005.
#  purpose:     Translation of lammps output into pdb format for use in
#               conjunction with vmd
#
#  Notes:       Copyright by author for Sandia National Laboratories
#    20040903   Conception date of v1.0: rudimentary script for collagen
#               project.
#    20050423   Conception date of v2.0:
#               - changed data access through indexing data directly on disk;
#               - added all command line options
#    20050425   Corrected focussing to use a target molecule's moment of
#               inertia to determine its principal axis: for computational
#               ease, a 3D orientational vector is calculated from two 2D
#               projections
#    20050701   Changed create_names() to loosen mass recognition criterion and
#               corrected created_atom_ids() to properly deal with the end of
#               the data stream.

# subroutines

  sub test
  {
    my $name            = shift(@_);

    printf("Error: file %s not found\n", $name) if (!scalar(stat($name)));
    return !scalar(stat($name));
  }


  sub initialize
  {
    my $k               = 0;
    my @options         = ("-help", "-nstart", "-dn", "-cut", "-repair",
                           "-units", "-quiet", "-nodetect", "-data", "-pbc",
                           "-focus", "-center", "-exclude");
    my @remarks         = ("display this message",
                           "starting frame [-1]",
                           "frames to skip when creating multiple frames [0]",
                           "cut bonds crossing over box edge [off]",
                           "repair bonds [off]",
                           "dump file entries have units [off]",
                           "turn display information off",
                           "turn auto-detection of masses off",
                           "use data file other than project.data",
                           "apply periodic boundary to molecules []",
                           "molecules to focus on []",
                           "molecules to use as center of mass []",
                           "exclude molecules from output []");
    my @notes           = (
      "* Multiple frames are processed when dn > 0.",
      "* Only the last frame is converted when nstart < 0.",
      "* Focussing superceedes centering.",
      "* Focussing uses the moment of inertia to determine the molecules'",
      "  principal axis; the molecules are rotated and translated to the lab",
      "  frame, using the focus molecules as reference.",
      "* Expected files in current directory: project.data, project.dump",
      "* Generated output files: project_trj.psf, project_trj.pdb",
      "* Uses project_ctrl.psf if available");


    $program            = "lammps2pdb";
    $version            = "2.2.5";
    $year               = "2007";
    $cut                = 0;
    $repair             = 0;
    $units              = 0;
    $info               = 1;
    $nstart             = -1;
    $dn                 = 0;
    $detect             = 1;

    foreach (@ARGV)
    {
      if (substr($_, 0, 1) eq "-")
      {
        my $k           = 0;
        my @arg         = split("=");
        my $switch      = (($arg[1] eq "")||($arg[1] eq "on")||($arg[1]!=0));
        $arg[0]         = lc($arg[0]);
        foreach (@options)
        {
          last if ($arg[0] eq substr($_, 0, length($arg[0])));
          ++$k;
        }
        $help           = 1 if (!$k--);
        $nstart         = $arg[1] if (!$k--);
        $dn             = $arg[1] if (!$k--);
        $cut            = $switch if (!$k--);
        $repair         = $switch if (!$k--);
        $units          = $switch if (!$k--);
        $info           = $switch ? 0 : 1 if (!$k--);
        $detect         = $switch ? 0 : 1 if (!$k--);
        $data_name      = $arg[1] if (!$k--);
        if (!$k--) {
          if ($switch) { $pbc{ALL} = 1; }
          else { foreach (split(",",$arg[1])) { $pbc{uc($_)} = 1; }}}
        if (!$k--) { foreach (split(",",$arg[1])) { $focus{uc($_)} = uc($_);}}
        if (!$k--) { foreach (split(",",$arg[1])) { $center{uc($_)} = uc($_);}}
        if (!$k--) { foreach (split(",",$arg[1])) { $exclude{uc($_)}=uc($_);}}
      }
      else { $project   = $_ if (!$k++); }
    }
    if (($k<1)||$help)
    {
      printf("%s v%s (c)%s by Pieter J. in \'t Veld for SNL\n\n",
        $program, $version, $year);
      printf("Usage:\n  %s.pl [-option[=#] ..] project\n\n", $program);
      printf("Options:\n");
      for (my $i=0; $i<scalar(@options); ++$i)
      {
        printf("  %-10.10s %s\n", $options[$i], $remarks[$i]);
      }
      printf("\nNotes:\n") if (scalar(@notes));
      foreach(@notes) { printf("  %s\n", $_); };
      printf("\n");
      exit(-1);
    }
    printf("%s v%s (c)%s\n\n", $program, $version, $year) if ($info);

    $data_name          = $project.".data" if ($data_name eq "");
    $traject_name       = $project.".dump";
    $pdb_name           = $project."_trj.pdb";
    $psf_name           = $project."_trj.psf";
    $psf_ctrl_name      = $project."_ctrl.psf";
    $psf_ctrl           = scalar(stat($psf_ctrl_name));

    $traject_flag       = scalar(stat($traject_name));
    my $flag            = test($data_name);
    printf("Conversion aborted\n\n") if ($flag);
    exit(-1) if ($flag);

    # data input

    create_atom_ids();
    create_bonds();
    open(TRAJECT, "<".$traject_name) if ($traject_flag);
    open(PSF_CTRL, "<".$psf_ctrl_name) if ($psf_ctrl);

    # open output

    open(PSF, ">".$psf_name) if (!$psf_ctrl);
    open(PDB, ">".$pdb_name);

    # align center with focus

    %center             = %focus if (scalar(%focus));
  }


# Vector routines

  sub V_String                                  # V_String(@a)
  {
    return "{".join(", ", @_)."}";
  }


  sub V_Round
  {
    foreach(@_) { $_ = ($_<0 ? -int(-$_*1e11+0.5) : int($_*1e11+0.5))/1e11; }
    return @_;
  }


  sub V_Add                                     # V_Add(@a, @b) = @a + @b;
  {
    return (@_[0]+@_[3], @_[1]+@_[4], @_[2]+@_[5]);
  }


  sub V_Subtr                                   # V_Subtr(@a, @b) = @a - @b;
  {
    return (@_[0]-@_[3], @_[1]-@_[4], @_[2]-@_[5]);
  }


  sub V_Dot                                     # V_Dot(@a, @b) = @a . @b;
  {
    return (@_[0]*@_[3]+@_[1]*@_[4]+@_[2]*@_[5]);
  }


  sub V_Cross                                   # V_Cross(@a, @b) = @a x @b;
  {
    return (@_[1]*@_[5]-@_[2]*@_[4], @_[2]*@_[3]-@_[0]*@_[5],
      @_[0]*@_[4]-@_[1]*@_[3]);
  }


  sub V_Mult                                    # V_Mult($f, @a) = $f * @a;
  {
    return (@_[0]*@_[1], @_[0]*@_[2], @_[0]*@_[3]);
  }


  sub V_Norm                                    # V_Norm(@a) = @a/|@a|;
  {
    return V_Mult(1/sqrt(V_Dot(@_[0,1,2],@_[0,1,2])), @_[0,1,2]);
  }


  sub pbc                                       # periodic -0.5*L <= x < 0.5*L
  {
    my $x               = @_[0]/@_[1]+0.5;

    return @_[0]-@_[1]*($x<0 ? int($x)-1 : int($x));
  }


  sub V_PBC                                     # V_PBC(@a, @l)
  {
    return (pbc(@_[0], @_[3]), pbc(@_[1], @_[4]), pbc(@_[2], @_[5]));
  }


  sub V_Cmp                                     # V_Cmp(abs(@a), abs(@b))
  {
    return -1 if (abs($_[0])<abs($_[3]));
    return  1 if (abs($_[0])>abs($_[3]));
    return -1 if (abs($_[1])<abs($_[4]));
    return  1 if (abs($_[1])>abs($_[4]));
    return -1 if (abs($_[2])<abs($_[5]));
    return  1 if (abs($_[2])>abs($_[5]));
    return 0;
  }


  sub V_Sort                                    # sort on descending absolute
  {                                             # value
    my @v               = @_;

    for (my $i=0; $i<scalar(@v)-3; $i+=3)
    {
      for (my $j=$i+3; $j<scalar(@v); $j+=3)
      {
        my @u           = @v[$i, $i+1, $i+2];
        next if (V_Cmp(@u, @v[$j, $j+1, $j+2])>=0);
        @v[$i,$i+1,$i+2]= @v[$j,$j+1,$j+2];
        @v[$j,$j+1,$j+2]= @u;
      }
    }
    return @v;
  }


# Matrix routines

  sub M_String                                  # M_String(@A)
  {
    my @M;

    for (my $i=0; $i<3; ++$i) { push(@M, V_String(splice(@_, 0, 3))); }
    return "{".join(", ", @M)."}";
  }


  sub M_Round
  {
    return V_Round(@_);
  }


  sub M_Transpose                               # M_Transpose(@A) = (@A)^t;
  {
    return (@_[0], @_[3], @_[6], @_[1], @_[4], @_[7], @_[2], @_[5], @_[8]);
  }


  sub M_Dot                                     # M_Dot(@A, @B) = @A . @B;
  {
    return (
      V_Dot(@_[0,1,2], @_[ 9,12,15]), V_Dot(@_[0,1,2], @_[10,13,16]),
        V_Dot(@_[0,1,2], @_[11,14,17]),
      V_Dot(@_[3,4,5], @_[ 9,12,15]), V_Dot(@_[3,4,5], @_[10,13,16]),
        V_Dot(@_[3,4,5], @_[11,14,17]),
      V_Dot(@_[6,7,8], @_[ 9,12,15]), V_Dot(@_[6,7,8], @_[10,13,16]),
        V_Dot(@_[6,7,8], @_[11,14,17]));
  }


  sub M_Det                                     # M_Det(@A) = | @A |
  {
    return V_Dot(@_[0,1,2], V_Cross(@_[3,4,5], @_[6,7,8]));
  }


  sub M_Mult                                    # M_Mult($a, @A) = $a * @A
  {
    return (
      @_[0]*@_[1], @_[0]*@_[2], @_[0]*@_[3],
      @_[0]*@_[4], @_[0]*@_[5], @_[0]*@_[6],
      @_[0]*@_[7], @_[0]*@_[8], @_[0]*@_[9]);
  }


  sub M_Unit { return (1,0,0, 0,1,0, 0,0,1); }

  sub PI { return 4*atan2(1,1); }

  sub M_Rotate                                  # M_Rotate($n, $alpha) = @R_$n;
  {                                             # vmd convention
    my $n               = shift(@_);
    my $alpha           = shift(@_)*PI()/180;
    my $cos             = cos($alpha);
    my $sin             = sin($alpha);

    $cos                = 0 if (abs($cos)<1e-16);
    $sin                = 0 if (abs($sin)<1e-16);
    return (1,0,0, 0,$cos,-$sin, 0,$sin,$cos) if ($n==0); # around x-axis
    return ($cos,0,$sin, 0,1,0, -$sin,0,$cos) if ($n==1); # around y-axis
    return ($cos,-$sin,0, $sin,$cos,0, 0,0,1) if ($n==2); # around z-axis
    return M_Unit();
  }


  sub M_RotateNormal                            # returns R.(1,0,0) = @a/|@a|;
  {
    my @R               = M_Unit();
    my @n               = V_Mult(1.0/sqrt(V_Dot(@_[0,1,2], @_)), @_);
    my $sina            = $n[1]<0 ? -sqrt($n[1]*$n[1]+$n[2]*$n[2]) :
                                     sqrt($n[1]*$n[1]+$n[2]*$n[2]);

    if ($sina)
    {
      my $cosa          = $n[0];
      my $cosb          = $n[1]/$sina;
      my $sinb          = $n[2]/$sina;
      @R                = (
        $cosa      ,-$sina      ,      0,
        $cosb*$sina, $cosb*$cosa, -$sinb,
        $sinb*$sina, $sinb*$cosa,  $cosb);
    }
    return @R;
  }


  sub MV_Dot                                    # MV_Dot(@A, @b) = @A . @b;
  {
    return (V_Dot(@_[0,1,2], @_[9,10,11]), V_Dot(@_[3,4,5], @_[9,10,11]),
      V_Dot(@_[6,7,8], @_[9,10,11]));
  }


  sub M_Sweep
  {
    my @M               = @_;

    for (my $i=0; $i<3; ++$i)
    {
      @M                = V_Sort(@M);
      my @v             = splice(@M, 3*$i, 3);
      if (abs($v[$i])>1e-10)
      {
        @v              = V_Mult(1/$v[$i], @v);
        @M[0,1,2]       = V_Subtr(@M[0,1,2], V_Mult($M[$i], @v));
        @M[3,4,5]       = V_Subtr(@M[3,4,5], V_Mult($M[$i+3], @v));
      }
      @M                = (@v, @M) if ($i==0);
      @M                = (@M[0,1,2], @v, @M[3,4,5]) if ($i==1);
      @M                = (@M, @v) if ($i==2);
    }
    return V_Round(@M);
  }


# Complex routines

  sub C_Add                                     # z = z1 + z2
  {
    return (@_[0]+@_[2], @_[1]+@_[3]);
  }


  sub C_Subtr                                   # z = z1 - z2
  {
    return (@_[0]-@_[2], @_[1]-@_[3]);
  }


  sub C_Mult                                    # z = z1 * z2
  {
    return (@_[0]*@_[2]-@_[1]*@_[3], @_[0]*@_[3]+@_[2]*@_[1]);
  }


  sub C_Div                                     # z = z1 / z2
  {
    my $num             = @_[2]*@_[2]+@_[3]*@_[3];

    return ((@_[0]*@_[2]+@_[1]*@_[3])/$num, (@_[1]*@_[2]-@_[0]*@_[3])/$num);
  }


  sub C_Pow                                     # z = z1 ^ a
  {
    my $r               = (sqrt(@_[0]*@_[0]+@_[1]*@_[1]))**@_[2];
    my $a               = @_[2]*atan2(@_[1], @_[0]);

    return ($r*cos($a), $r*sin($a));
  }


  sub C_Correct
  {
    return (abs(@_[0])<1e-14 ? 0 : @_[0], abs(@_[1])<1e-14 ? 0 : @_[1]);
  }


  sub C_Zero
  {
    return (abs(@_[0])<1e-8 ? 0 : @_[0], abs(@_[1])<1e-8 ? 0 : @_[1]);
  }


  sub C_Conj
  {
    return (@_[0], -@_[1]);
  }


  sub C_String
  {
    return @_[0]." + ".@_[1]."i";
  }


# analytical roots

  sub R_Sort
  {
    my $n               = scalar(@_);

    for (my $i=0; $i<$n-2; $i+=2)
    {
      for (my $j=$i+2; $j<$n; $j+=2)
      {
        if (@_[$j]<@_[$i]) {
          my @t = @_[$i,$i+1]; @_[$i,$i+1] = @_[$j,$j+1]; @_[$j,$j+1] = @t; }
        else { if ((@_[$j]==@_[$i])&&(@_[$j+1]<@_[$i+1])) {
          my @t = @_[$i,$i+1]; @_[$i,$i+1] = @_[$j,$j+1]; @_[$j,$j+1] = @t; } }
      }
    }
    return @_;
  }


  sub R_First
  {
    return (0, 0) if (abs(@_[1])<1e-14);
    return (-@_[0]/@_[1], 0);
  }


  sub R_Second
  {
    return R_First(@_) if (abs(@_[2]<1e-14));
    my @z               = (-@_[1]/@_[2]/2.0, 0);
    my @root            = C_Correct(C_Pow(($z[0]*$z[0]-@_[0]/@_[2], 0), 1/2));
    return (C_Zero(C_Subtr(@z, @root)), C_Zero(C_Add(@z, @root)));
  }


  sub R_Third
  {
    return R_Second(@_) if (abs(@_[3])<1e-14);
    my $c3              = 3*@_[3];
    my $f1              = @_[1]*$c3-@_[2]*@_[2];
    my $f2              = 0.5*@_[2]*(3*$f1+@_[2]*@_[2])-1.5*@_[0]*$c3*$c3;
    my $f3              = -@_[2]/$c3;
    my @A               = (0, 0);
    my @B               = (0, 0);
    my @z1              = (0.5, 0.5*sqrt(3));
    my @z2              = C_Conj(@z1);

    if (abs($f1)<1e-3) {                                # limit f1 -> 0
      @A                = ($f2<0 ? abs(2*$f2)**(1/3) : 0, 0); }
    else {
      if (abs($f2)<1e-14) {                             # limit f2 -> 0
        my $f           = sqrt(abs($f1))/$c3;
        @A              = $f1<0 ? (-$f*$z1[1], 0.5*$f) : ($f, 0);
        @B              = $f1<0 ? (-$A[0], $A[1]) : ($f, 0); }
      else {
        @B              = C_Pow(C_Add(($f2, 0),
                            C_Pow(($f1*$f1*$f1+$f2*$f2, 0), 1/2)), 1/3);
        @A              = C_Div(($f1/$c3, 0), @B);
        @B              = ($B[0]/$c3, $B[1]/$c3); } }
    return R_Sort(
      C_Zero(C_Add(($f3, 0), C_Subtr(C_Mult(@z1, @A), C_Mult(@z2, @B)))),
      C_Zero(C_Add(($f3, 0), C_Subtr(C_Mult(@z2, @A), C_Mult(@z1, @B)))),
      C_Zero(C_Add(($f3, 0), C_Subtr(@B, @A))));
  }


# general read file operators

  sub advance_read
  {
    my $input           = shift;
    my $dlines          = shift;
    my $read;

    while ($dlines--) { $read = <$input>; }
    return $read;
  }


  sub rewind_read
  {
    my $input           = shift;

    seek($input, 0, SEEK_SET);
  }


# create crossreference tables

  sub create_psf_index                                  # make an index of id
  {                                                     # locations
    my @psf_ids         = ("!NATOM","!NBOND:","!NTHETA:","!NPHI:","!NIMPHI:");
    my @ids             = (atoms, bonds, angles, dihedrals, impropers);
    my $k               = 0;
    my %hash;

    printf("Info: creating PSF index\n") if ($info);
    foreach (@psf_ids) { $hash{$_} = shift(@ids); };
    seek(PSF_CTRL, 0, SEEK_SET);
    while (<PSF_CTRL>)
    {
      chop();
      my @tmp           = split(" ");
      my $n             = $hash{$tmp[1]};
      $PSFIndex{$n}     = tell(PSF_CTRL)." ".$tmp[0] if ($n ne "");
    }
  }


  sub psf_goto                                          # goto $ident in <PSF>
  {
    create_psf_index() if (!scalar(%PSFIndex));
    my $id              = shift(@_);
    my @n               = split(" ", $PSFIndex{$id});

    @PSFBuffer          = ();
    if (!scalar(@n))
    {
      printf("Warning: PSF index for $id not found\n");
      seek(PSF_CTRL, 0, SEEK_END);
      return -1;
    }
    seek(PSF_CTRL, $n[0], SEEK_SET);
    return $n[1];
  }


  sub psf_get
  {
    @PSFBuffer          = split(" ", <PSF_CTRL>) if (!scalar(@PSFBuffer));
    return splice(@PSFBuffer, 0, shift(@_));
  }


  sub create_data_index
  {
    my $init            = 3;
    my @tmp;
    my $id;
    my %hash;
    my %size;

    foreach ((masses,atoms,bonds,angles,dihedrals,impropers)) { $hash{$_}=$_; }
    open(DATA, "<".$data_name);
    my $tmp = <DATA>;                                   # skip first line
    while ($init&&!eof(DATA))                           # interpret header
    {
      @tmp              = split(" ", <DATA>);
      --$init if (!scalar(@tmp));
      foreach (@tmp) { $_ = lc($_); }
      $tmp[1]           = masses if ($tmp[1]." ".$tmp[2] eq "atom types");
      $L[0]             = $tmp[0]." ".($tmp[1]-$tmp[0])
                                        if (join(" ",@tmp[2,3]) eq "xlo xhi");
      $L[1]             = $tmp[0]." ".($tmp[1]-$tmp[0])
                                        if (join(" ",@tmp[2,3]) eq "ylo yhi");
      $L[2]             = $tmp[0]." ".($tmp[1]-$tmp[0])
                                        if (join(" ",@tmp[2,3]) eq "zlo zhi");
      if (($id = $hash{$tmp[1]}) ne "") { $size{$id}    = $tmp[0]; }
    }
    @l                  = ();
    for (my $i=0; $i<3; ++$i)
    {
      @tmp              = split(" ", $L[$i]);
      $l[$i]            = $tmp[0];
      $l[$i+3]          = $tmp[1];
    }
    while (!eof(DATA))                                  # interpret data
    {
      @tmp              = split(" ", <DATA>);
      my $skip          = <DATA> if (($id = $hash{lc($tmp[0])}) ne "");
      $DATAIndex{$id}   = tell(DATA)." ".$size{$id} if ($id ne "");
    }
  }


  sub goto_data
  {
    create_data_index() if (!scalar(%DATAIndex));
    my $id              = shift(@_);
    my @n               = split(" ", $DATAIndex{$id});

    if (!scalar(@n))
    {
      printf("Warning: DATA index for $id not found\n");
      seek(DATA, 0, SEEK_END);
      return -1;
    }
    seek(DATA, $n[0], SEEK_SET);
    return $n[1];
  }


# create atom and residue identifiers

  sub create_names
  {
    return if (scalar(@names));
    my $n               = goto_data(masses);
    my @mass            = (1, 12, 14, 16, 32.1, 4, 20.2, 40.1, 65.4,
                           55.8, 31, 35.5, 23, 24.3, 39.1, 28.1);
    my @name            = ("H", "C", "N", "O", "S", "HE", "NE", "CA",
                           "ZN", "FE", "P", "CL", "NA", "MG", "K", "SI");
    my @unknown         = ("X", "Y", "Z");
    my @letter          = ("X", "Y", "Z", "P", "Q", "R", "S", "A", "B", "C");
    my $k               = 0;
    my %atom;

    $names[0]           = "";
    foreach (@mass) { $atom{$_} = shift(@name); }
    for (my $i=1; $i<=$n; ++$i)
    {
      my @tmp           = split(" ", <DATA>);
      my $j             = $tmp[0];
      $tmp[1]           = int($tmp[1]*10+0.5)/10;
      if ((($names[$j] = $atom{$masses[$j] = $tmp[1]}) eq "")||!$detect)
      {
        $names[$j]      = $letter[$k].$unknown[0];
        if (++$k>9) { $k = 0; shift(@unknown); }
      }
    }
  }


  sub create_position
  {
    my @data            = @_;
    my $p               = $data[1]." ".$data[2];
    my $k;

    for ($k=0; ($k<scalar(@data))&&(substr($data[$k],0,1) ne "#"); ++$k) { }
    @data               = splice(@data, $k-($k<8 ? 3 : 6), $k<8 ? 3 : 6);
    foreach (@L)
    {
      my @l             = split(" ");
      $p                = join(" ", $p, ($data[0]-$l[0])/$l[-1]+$data[3]);
      shift(@data);
    }
    return $p;
  }


  sub create_atom_ids
  {
    my $res             = 0;
    my $last            = 1;
    my @data;
    my $n;
    my $k;
    my $tmp;
    my $id;
    my %link;
    my %special;

    printf("Info: creating atom ids\n") if ($info);
    create_names();
    $n                  = goto_data(atoms);
    foreach ((CL, NA, MG, K)) { $special{$_} = "ION"; }
    $special{"H H O"}   = "HOH";
    for (my $i=0; $i<=$n; ++$i)
    {
      @data             = split(" ", <DATA>) if ($i<$n);
      $traject[$i]      = create_position(@data) if ($i<$n);    # positions
      if ($i<$n ? ($mol[$i+1] = $data[1])!=$last : 1)           # residues
      {
        if ((($tmp = $link{$id = join(" ", sort(split(" ", $id)))}) eq "")&&
            (($tmp = $special{$id}) eq ""))
        {
          $link{$id}    =
              $tmp      = "R".($res<10 ? "0" : "").$res;
          ++$res;
        }
        $cluster[$last] = "PROT";
        $cluster[$last] = "WATR" if ($tmp eq "HOH");
        $cluster[$last] = "SALT" if ($tmp eq "ION");
        $residue[$last] = $tmp;
        $last           = $data[1];
        $id             = "";
      }
      $id               = join(" ", $id, $names[$data[2]]);
    }
  }


  sub crossover
  {
    my @d               = V_Subtr((split(" ", $position[@_[0]]))[0,1,2],
                                  (split(" ", $position[@_[1]]))[0,1,2]);

    $d[0]               /= $l[3];
    $d[1]               /= $l[4];
    $d[2]               /= $l[5];
    return (($d[0]<-0.5)||($d[0]>=0.5)||($d[1]<-0.5)||($d[1]>=0.5)||
      ($d[2]<-0.5)||($d[2]>=0.5));
  }


  sub delete_crossovers
  {
    my $n               = scalar(@bonds);
    my $i               = 0;

    printf("Info: deleting crossover bonds\n") if ($info);
    while ($i<$n)                                       # take out crossovers
    {
      if (crossover(split(" ", $bonds[$i]))) { splice(@bonds, $i, 1); --$n; }
      else { ++$i; }
    }
  }


  sub delete_exclude
  {
    my $n               = scalar(@bonds);
    my $i               = 0;

    printf("Info: deleting excluded bonds\n") if ($info);
    while ($i<$n)
    {
      my $m             = $mol[(split(" ", $bonds[$i]))[0]+1];
      if (($exclude{$m} ne "")||($exclude{$residue[$m]} ne "")
        ||($exclude{$cluster[$m]} ne "")) { splice(@bonds, $i, 1); --$n; }
      else { ++$i; }
    }
  }


  sub create_bonds
  {
    my $n               = goto_data(bonds);

    printf("Info: creating bonds\n") if ($info);
    for (my $i=0; $i<$n; ++$i)
    {
      my @arg           = split(" ", <DATA>);
      $bonds[$i]        = ($arg[2]-1)." ".($arg[3]-1);
    }
  }


# traject operators

  sub advance_traject
  {
    my $subject         = "item: ".lc(shift(@_));
    my $advance         = 1;

    while (!eof(TRAJECT)&&(substr(lc(join(" ", split(" ", <TRAJECT>))),
                                  0,length($subject)) ne $subject)) {}
  }


  sub read_traject
  {
    my @box;
    my @l;

    advance_traject("timestep");
    my $timestep        = (split(" ", <TRAJECT>))[0];
    advance_traject("number of atoms");
    my $n               = (split(" ", <TRAJECT>))[0];
    advance_traject("box bounds");
    for (my $i=0; $i<3; ++$i)
    {
      my @data          = split(" ", <TRAJECT>);
      $box[$i]          = $data[0];                     # box edge
      $l[$i]            = $data[1]-$data[0];            # box length
    }
    advance_traject("atoms");
    for (my $i=0; $i<$n; ++$i)
    {
      my @data          = split(" ", <TRAJECT>);        # read data
      $traject[$data[0]-1] = join(" ", @data);          # store data in order
    }
    return ($timestep, $n, @box, @l);
  }


# pdb format

  sub eigen_vector                              # eigen_vector(@A, $l)
  {
    my @A               = splice(@_,0,9);
    my $l               = shift(@_);

    $A[0]               -= $l;
    $A[4]               -= $l;
    $A[8]               -= $l;
    @A                  = M_Sweep(@A);
    return V_Norm(-$A[2], -$A[5], 1) if (($A[0]==1)&&($A[4]==1));
    return V_Norm(-$A[1], 1, -$A[7]) if (($A[0]==1)&&($A[8]==1));
    return V_Norm(1, -$A[3], -$A[6]) if (($A[4]==1)&&($A[8]==1));
    return (1,0,0) if ($A[0]==1);
    return (0,1,0) if ($A[4]==1);
    return (0,0,1) if ($A[8]==1);
    return (0,0,0);
  }


  sub pdb_inertia
  {
    my @s               = (
      @_[3]-@_[0]*@_[0], @_[4]-@_[1]*@_[1], @_[5]-@_[2]*@_[2],
      @_[6]-@_[0]*@_[1], @_[7]-@_[0]*@_[2], @_[8]-@_[1]*@_[2]);
    my @c               = (
      ($s[0]+$s[1])*(($s[0]+$s[2])*($s[1]+$s[2])-$s[3]*$s[3])           # c0
        -$s[5]*($s[3]*$s[4]+($s[1]+$s[2])*$s[5])
        -$s[4]*(($s[0]+$s[2])*$s[4]+$s[3]*$s[5]),
      -($s[0]+$s[2])*($s[1]+$s[2])-($s[0]+$s[1])*($s[0]+$s[1]+2*$s[2])  # c1
        +$s[3]*$s[3]+$s[4]*$s[4]+$s[5]*$s[5],
      2*($s[0]+$s[1]+$s[2]),                                            # c2
      -1);                                                              # c3
    my @sol             = R_Third(@c);
    my @M               = ($s[1]+$s[2], -$s[3], -$s[4],
                           -$s[3], $s[0]+$s[2], -$s[5],
                           -$s[4], -$s[5], $s[0]+$s[1]);
    my @a               = eigen_vector(@M, $sol[0]);
    my @b               = eigen_vector(@M, $sol[2]);
    my @A               = M_Transpose(M_RotateNormal(@a));
    my @B               = M_Dot(M_RotateNormal(0,1,0),
                            M_Transpose(M_RotateNormal(MV_Dot(@A,@b))));
    return M_Dot(@B, @A);
  }



  sub pdb_focus         # using moment of inertia
  {
    my @R               = pdb_inertia(@_);

    printf("Info: focussing\n") if ($info);
    foreach (@position)
    {
      my @p             = split(" ");
      $_                = join(" ", MV_Dot(@R, @p[0,1,2]), @p[3,4]);
    }
  }


  sub pdb_center        
  {
    my @c               = splice(@_, 0, 3);

    printf("Info: centering\n") if ($info);
    @l[0,1,2]           = V_Mult(-1/2, @l[3,4,5]);
    foreach (@position)
    {
      my @p             = split(" ");
      $_                = join(" ", V_Subtr(@p[0,1,2], @c), @p[3,4]);
    }
  }


  sub pdb_pbc
  {
    printf("Info: applying periodicity\n") if ($info);
    foreach (@position)
    {
      my @p             = split(" ");
      my $m             = $mol[$p[3]];
      $_                = join(" ", V_PBC(@p[0,1,2], @l[3,4,5]), @p[3,4])
            if ($pbc{ALL}||$pbc{$m}||$pbc{$residue[$m]}||$pbc{$cluster[$m]});
    }
  }


  sub pdb_repair_bonds
  {
    return if (!$pbc);
    printf("Info: repairing bonds\n") if ($info);
    foreach (@bonds)
    {
      my @b             = split(" ");
      my @p             = split(" ", $position[$b[0]]);
      my @q             = split(" ", $position[$b[1]]);
      $position[$b[1]]  = join(" ", V_Add(@p[0,1,2], V_PBC(
                          V_Subtr(@q[0,1,2], @p[0,1,2]), @l[3,4,5])), @q[3,4]);
    }
  }


  sub pdb_positions
  {
    my @m               = (0,0,0,0,0,0,0,0,0);
    my $nmass           = 0;
    my $i               = 0;
    my $mass;
    my @p;
    my $d;

    foreach (@traject)
    {
      my @arg           = split(" ");
      my $m             = $mol[$arg[0]];
      next if (($exclude{$m} ne "")||($exclude{$residue[$m]} ne "")
        ||($exclude{$cluster[$m]} ne ""));
      if ($units)
      {
        $p[0]           = $arg[2]+$arg[5]*$l[3];
        $p[1]           = $arg[3]+$arg[6]*$l[4];
        $p[2]           = $arg[4]+$arg[7]*$l[5];
      }
      else
      {
        $p[0]           = $l[0]+($arg[2]+$arg[5])*$l[3];
        $p[1]           = $l[1]+($arg[3]+$arg[6])*$l[4];
        $p[2]           = $l[2]+($arg[4]+$arg[7])*$l[5];
      }
      $position[$i++]   = join(" ", @p, $arg[0], $arg[1]);
      next if (($center{$m} eq "")&&($center{$cluster[$m]} eq ""));
      $nmass            += $mass = $masses[$arg[1]];    # inertia necessities:
      $m[0]             += $d = $mass*$p[0];            # <x>
      $m[3]             += $d*$p[0];                    # <x^2>
      $m[6]             += $d*$p[1];                    # <xy>
      $m[7]             += $d*$p[2];                    # <xz>
      $m[1]             += $d = $mass*$p[1];            # <y>
      $m[4]             += $d*$p[1];                    # <y^2>
      $m[8]             += $d*$p[2];                    # <yz>
      $m[2]             += $d = $mass*$p[2];            # <z>
      $m[5]             += $d*$p[2];                    # <z^2>
    }
    pdb_center(M_Mult(1/$nmass, @m)) if ($nmass);
    pdb_focus(M_Mult(1/$nmass, @m)) if ($nmass && scalar(%focus));
    pdb_pbc() if (scalar(%pbc));
    pdb_repair_bonds() if ($repair);
  }


  sub pdb_header
  {
    printf(PDB "REMARK  \n");
    printf(PDB "REMARK  ".$project."_trj.pdb GENERATED FROM ".$project.
      ($psf_ctrl ? "_ctrl.psf" : ".data")." AND $project.dump\n");
    printf(PDB "REMARK  CREATED BY $program v$version ON %s", `date`);
    printf(PDB "REMARK  \n");
    printf(PDB "CRYST1 ");
    printf(PDB "%8.8s ", $l[3]);
    printf(PDB "%8.8s ", $l[4]);
    printf(PDB "%8.8s ", $l[5]);
    printf(PDB "%6.6s ", "90");
    printf(PDB "%6.6s ", "90");
    printf(PDB "%6.6s ", "90");
    printf(PDB "%-11.11s", "P 1");
    printf(PDB "%3.3s\n", "1");
  }


  sub pdb_atoms
  {
    my $n               = 0;
    my $l               = 0;
    my $last            = 0;
    my @base;

    pdb_positions();
    psf_goto(atoms) if ($psf_ctrl);
    printf("Info: writing positions for timestep %d\n", $description[0])
      if ($info);
    foreach (@position)
    {
      my @p             = split(" ");
      my $nres          = $mol[$p[3]];
      my @psf           = split(" ", advance_read(PSF_CTRL, $p[3]-$l))
                                                              if ($psf_ctrl);
      @base             = @p[0,1,2] if ($last!=$nres);
      #@p               = V_Add(V_PBC(V_Subtr(@p, @base), @l[3,4,5]), @l);
      foreach (@p) { $_ = 0 if (abs($_)<1e-4); }
      printf(PDB "ATOM  ");                             # pdb command
      printf(PDB "%6.6s ", ++$n);                       # atom number
      printf(PDB "%-3.3s ",
        $psf_ctrl ? $psf[4] : $names[$p[4]]);           # atom name
      printf(PDB "%-3.3s ",
        $psf_ctrl ? $psf[3] : $residue[$nres]);         # residue name
      printf(PDB "%5.5s ", $nres);                      # residue number
      printf(PDB "%3.3s ", "");                         # empty placeholder
      printf(PDB "%7.7s %7.7s %7.7s ",
        $p[0], $p[1], $p[2]);                           # positions
      printf(PDB "%5.5s %5.5s %4.4s ",
        "1.00", "0.00", "");                            # trailing variables
      printf(PDB "%-4.4s\n",
        $psf_ctrl ? $psf[1] : $cluster[$nres]);         # cluster name
      $last             = $nres;
      $l                = $p[3];
    };
    printf(PDB "END\n");
  }


  sub pdb_timestep
  {
    pdb_header();
    pdb_atoms();
    return ();
  }


# psf format

  sub psf_header
  {
    printf(PSF "PSF\n");
    printf(PSF "\n");
    printf(PSF "%8.8s !NTITLE\n", 2);
    printf(PSF "REMARK   ".$project."_trj.psf GENERATED FROM $project.data\n");
    printf(PSF "REMARK   CREATED BY $program v$version ON %s", `date`);
    printf(PSF "\n");
  }


  sub psf_atoms
  {
    my $n               = goto_data(atoms);
    my $natom           = 0;
    my $l               = 0;
    my $k               = 0;
    my @extra;

    for (my $i=0; $i<$n; ++$i)
    {
      my @arg           = split(" ", <DATA>);
      if (!$k) {
        for ($k=0; ($k<scalar(@arg))&&(substr($arg[$k],0,1) ne "#"); ++$k) {} }
      $extra[$arg[0]-1] = $arg[1]." ".$arg[2]." ".$arg[3];
    }
    printf(PSF "%8.8s !NATOM\n", scalar(@position));
    foreach (@position)
    {
      my @p             = split(" ");
      my @arg           = split(" ", $extra[$natom]);
      printf(PSF "%8.8s ", ++$natom);                   # atom number
      printf(PSF "%-4.4s ", $cluster[$arg[0]]);         # cluster name
      printf(PSF "%-4.4s ", $arg[0]);                   # residue number
      printf(PSF "%-4.4s ", $residue[$arg[0]]);         # residue name
      printf(PSF "%-4.4s ", $names[$p[4]]);             # atom name
      printf(PSF "%4.4s ", $arg[1]);                    # atom number
      printf(PSF "%10.10s ", $k%3 ? $arg[2] : 0);       # atom charge
      printf(PSF "%4.4s ", "");                         # blank entry
      printf(PSF "%8.8s ", $masses[$arg[1]]);           # atom mass
      printf(PSF "%11.11s\n", "0");                     # trailing variable
      $l                = $p[3];
      last if ($natom>=$n)
    }
    printf(PSF "\n");
  }


  sub psf_bonds
  {
    my $npairs          = 0;

    delete_exclude() if (scalar(%exclude)>0);
    delete_crossovers() if ($cut);
    printf(PSF "%8.8s !NBOND\n", scalar(@bonds));
    foreach (@bonds)
    {
      my @arg           = split(" ");
      printf(PSF " %7.7s %7.7s", $arg[0]+1, $arg[1]+1);
      if (++$npairs>=4) { printf(PSF "\n"); $npairs = 0; }
    }
    if ($npairs>0) { printf(PSF "\n"); }
    printf(PSF "\n");
  }


# main

  initialize();

  # create .pdb file

  $ncurrent             = -1;
  while ($traject_flag&&!eof(TRAJECT))
  {
    @description        = read_traject();
    @l                  = splice(@description, -6, 6);
    next if (($nstart<0)||(++$ncurrent<$nstart));
    if ($dn<1) { pdb_timestep(); last; }
    $ncurrent           = pdb_timestep() if ($nstart||!($ncurrent%$dn));
    $nstart             = 0;
  }
  pdb_timestep() if ($nstart<0);

  # create .psf file

  if (!$psf_ctrl)
  {
    psf_header();
    psf_atoms();
    psf_bonds();
  }
  else
  {
    system("rm -f ".$psf_name);
    system("ln -s ".$psf_ctrl_name." ".$psf_name);
  }

  # add tail to files

  #printf(PDB "END\n");
  printf("\n") if ($info);
  close(PDB);
  close(PSF);

  close(TRAJECT);
  close(DATA);