Hi,
On Wed, Jun 21, 2023 at 5:29 PM Eduard Neu via AMBER <amber.ambermd.org> wrote:
> Basically, I want to determine every frame in which an oxygen atom of a
> water molecule is within 1 Ã… of the original crystal water.
It's a very interesting problem, and one that might be worth creating
a dedicated Action for. For now, I think I've come up with a somewhat
hacky solution that uses the 'grid' action to count waters in specific
positions. Caveat: I haven't really tested this in production so try
it out on a few trajectory frames first to see if your results make
sense.
I think you can do it in three passes and one final post-processing step:
1) Use the 'bounds' action to determine the size grid you need for
subsequent steps.
2) Use 'grid' on the reference frame to get a 'mask' you can use to
highlight the grid bins you want (corresponding to the
crystallographic waters).
3) Use 'grid' on the target trajectory to count waters in each bin,
making sure you rms fit to the reference structure.
4) Multiply the reference grid by the target grid to filter only the
positions you want.
The input would look something like this (note I'm using :HOH,WAT
since presumably in a PDB your waters are named HOH and in an Amber
topology they are named WAT):
# Load topology
parm 2b5t.pdb
# This mask will be used to rms fit
set FITMASK = ^1-4&.CA
# Replace these with your reference frame and target traj names
set REFTRAJ = 2b5t.pdb
set TGTTRAJ = 2b5t.pdb
# Load reference frame as a reference for rms
reference $REFTRAJ
# Create the common grid, center on origin
trajin $REFTRAJ
trajin $TGTTRAJ
autoimage origin
rms reference $FITMASK
bounds :HOH,WAT name MyGrid dx 1.0 out bounds.dat
run
writedata mygrid.dx MyGrid
clear trajin
# Determine water positions in Reference frame
readdata mygrid.dx name RefGrid
trajin $REFTRAJ
autoimage origin
rms reference $FITMASK
grid :HOH,WAT.O data RefGrid out refgrid.dx
run
clear trajin
# Count waters in target Trajectory
readdata mygrid.dx name TgtGrid
trajin $TGTTRAJ
autoimage origin
rms reference $FITMASK
grid :HOH,WAT.O data TgtGrid out tgtgrid.dx
run
clear trajin
# Multiply reference grid by target to only get positions in reference
Result = RefGrid * TgtGrid
# Write the results to a 'dat' file with 'sparse' keyword to ignore
zero occupancies
writedata Result.dat Result sparse
quit
The final results in Result.dat should have the XYZ position of each
relevant bin along with the non-normalized count in each bin (which
should just be the number of frames that bin was occupied). This is a
pretty complicated script but I think it should work - however, note
that whatever you use for FITMASK ** has to select the same number of
atoms in REFTRAJ and TGTTRAJ **. If your PDB is missing residues or
something this may not work. Anyway, give it a try and let me know if
you run into problems. Good luck!
-Dan
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jun 22 2023 - 07:30:02 PDT
