lammps2pdb.pl
00aef0fe0014lammm-devel

lammps2pdb.pl
View Options

			#!/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.5L <= x < 0.5L
			{
			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 ($rcos($a), $rsin($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);

lammps/tools/ch2lmp/lammps2pdb.pl00aef0fe0014lammm-devel

lammps2pdb.plView Options

lammps/tools/ch2lmp/lammps2pdb.pl
00aef0fe0014lammm-devel

lammps2pdb.pl
View Options