Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73374196
lammps2pdb.pl
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jul 20, 22:58
Size
27 KB
Mime Type
text/x-perl
Expires
Mon, Jul 22, 22:58 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18290084
Attached To
rLAMMPS lammps
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
3
D
orientational
vector
is
calculated
from
two
2
D
#
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
(
-
$
_
*
1
e11
+
0
.
5
)
:
int
(
$
_
*
1
e11
+
0
.
5
))
/
1
e11
;
}
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
)
<
1
e
-
16
);
$
sin
=
0
if
(
abs
(
$
sin
)
<
1
e
-
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
])
>
1
e
-
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
])
<
1
e
-
14
?
0
:
@
_
[
0
],
abs
(
@
_
[
1
])
<
1
e
-
14
?
0
:
@
_
[
1
]);
}
sub
C_Zero
{
return
(
abs
(
@
_
[
0
])
<
1
e
-
8
?
0
:
@
_
[
0
],
abs
(
@
_
[
1
])
<
1
e
-
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
])
<
1
e
-
14
);
return
(
-@
_
[
0
]
/@
_
[
1
],
0
);
}
sub
R_Second
{
return
R_First
(
@
_
)
if
(
abs
(
@
_
[
2
]
<
1
e
-
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
])
<
1
e
-
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
)
<
1
e
-
3
)
{
#
limit
f1
->
0
@
A
=
(
$
f2
<
0
?
abs
(
2
*
$
f2
)
**
(
1
/
3
)
:
0
,
0
);
}
else
{
if
(
abs
(
$
f2
)
<
1
e
-
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
(
$
_
)
<
1
e
-
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
);
Event Timeline
Log In to Comment