forked from lijiext/lammps
1138 lines
28 KiB
Raku
1138 lines
28 KiB
Raku
#!/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);
|
|
for (my $i=0; $i<2; ++$i) { my $tmp = <DATA>; } # skip first two lines
|
|
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);
|
|
|