Hi
On Mon, Jun 26, 2023 at 1:48 PM Eduard Neu <eduard.neu.fau.de> wrote:
> I was a little confused because I had 20,000 frames and the crystallographic waters remained quite stable during the MD but the results do not indicate that. I thought the results would be around 15,000. Also, the output included only 5 bins while I thought there should be at least 18. But maybe some of them are just empty. Additionally, the coordinates in the results file are not exactly the coordinates of the crystallographic waters.
So there will be some sensitivity here due to artifacts introduced by
grid bin size, i.e. if the water is just outside a voxel, or on the
border of a voxel, it will get divided between those bins. You could
try playing around with a slightly larger grid size to see if the
results look better. Another potential way to ameliorate this would be
to use the 'volmap' command instead of 'grid', and use the 'calcpeaks'
keyword to pick out regions of high water density, but probably this
would require a dedicated action to do "correctly".I'm travelling
right now, but it's something I can look into in the future .
-Dan
>
> I played around with the script and changed the MASK and the distance for the grid generation, the MASK for the reference grid and the MASK for the water count. I could generate results that were more in agreement with the trajectory but I was not sure if that is the right way to go.
>
> Do you have any idea how the script should be changed to generate more accurate results?
>
> Best and thanks again
> Eddy
>
> On 6/22/2023 4:00 PM, Daniel Roe wrote:
>
> 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 29 2023 - 06:30:03 PDT
