Dan, thanks for the effort and time given. The problem was resolved by carefully making topology, matching with the new coordinate.
--- On Mon, 23/8/10, Daniel Roe <daniel.r.roe.gmail.com> wrote:
From: Daniel Roe <daniel.r.roe.gmail.com>
Subject: Re: [AMBER] ptraj "closest"
To: "AMBER Mailing List" <amber.ambermd.org>
Date: Monday, 23 August, 2010, 10:07 PM
Hi,
What version of ptraj are you using and are all bugfixes applied?
On Mon, Aug 23, 2010 at 7:10 AM, Amor San Juan <amorsanjuan.yahoo.com> wrote:
> I wanted to select water molecules near the ligand, and then use this as input for mmpbsa calculation. Upon using the closest command of ptraj (refer to close.in below), it gave me an output with 10 selected water molecules at nearest distance to ligand. The problem was that the output trajectory contains all the tip3p waters, and I tried to strip this using ptraj1.in and then, use ptraj2.in to get new prmtop. Whenever I viewed in vmd, there are some long lines appearing.
My tests show that the 'closest' command is working as expected in AT
v1.4 ptraj. First you state that ptraj gave you "an output with 10
selected water molecules at nearest distance to ligand" (which seems
OK), but then you say "the output trajectory contains all the tip3p
waters". Which one is it? According to your 'close1.in' file, your
out1.mdcrd should be an Amber Trajectory that contains the coordinates
of the solute and 10 waters. The easiest way to check that this is
actually happening is to write a test PDB file:
trajin in.mdcrd
trajout test.pdb pdb
closest 10 :200
Assuming all is well so far, you can use your 'close1.in' input file
with ptraj to generate an Amber Trajectory. However, keep in mind that
the trajectory you generate has a different topology than your input
trajectory. Therefore you will have to generate a new topology file
that corresponds to the output trajectory (i.e. it has parameters for
the solute and 10 water molecules only) before you can read it back
into ptraj (otherwise you will get errors like > Mask [:211-9199]
represents 0 atoms !!!NO ATOMS DETECTED!!!). This can be done by
loading the test pdb you just generated containing the target number
of waters into tleap and using the 'set UNIT box' command to fix the
box info (see the AmberTools 1.4 manual section 3.4.38), then saving
the new topology with 'saveamberparm'.
Hopefully this information helps. Let me know if it is unclear or if
you have any more issues.
-Dan
>
> For reference, the :200 represents the ligand (with phosphate gp), :1-199 are residues and 210-9199 are tip3p water.
>
> A) Strategy-1
> close1.in
> -----
> trajin in.mdcrd
> trajout out1.mdcrd
> closest 10 :200
> image :WAT byres
> rms first mass :1-199
> -----
>
> ptraj1.in
> ----
> trajin out1.mdcrd
> strip :211-9199
> center :1-210
> image center
> trajout strip.mdcrd
> ----
>
> ptraj2.in
> ---
> trajin out1.mdcrd 1 1 1
> trajout top PDB
> center :1-210
> strip :211-9199
> ---
>
> The closest command in ptraj has the logical effect that it modifies the state of the system, ie., new numbering of atoms. Although in my output 'out.mdcrd' the 10 closest waters appeared immediately next to the solute, the remaining waters (211-9199) follows it. I do not want this remaining waters to be included for mmpbsa calculations. The attempt I did above is a strategy of first implementing the closest command, then the next ptraj script to strip the unwanted waters (far distant from ligand). Another strategy I did was to issue both closest and strip in one single script as follows:
>
> B)Strategy-2
>
> close2.in
> -------
> trajin in.mdcrd
>
> trajout out1.mdcrd
>
> closest 10 :200
>
> image :WAT byres
>
> rms first mass :1-199
> strip :211-9199
> -----
>
> Here is the error:
>
> Understanding how closest command works, it does change the state of the system, particularly the numbering, as evidenced of its inability to detect the mask 211-9199. If that is the case, how to strip the unwanted waters?
>
> My assumption is that since the closest command designates new atom numbering, using the strip command will not detect any specified water atoms that I desired to remove from the trajectory. This could also be the reason why am I viewing unwanted long lines in vmd.
>
> Alternatively, I read from the forum about this simple script:
>
> C)Strategy-3
> ----
> trajin out1.mdcrd
> strip "(:200 >.5.0)"
>
> trajout strip.mdcrd
> -----
>
> The script strip all atoms further than 5 Angstroms from
> the ligand (:200). Problem with this is that my input has the protein + ligand, so what the script does is it strip all atoms (including protein and water) located further than 5 angstroms.
>
> On the lighter side, apart from those 3 strategies I was able to view in vmd the 'out1.mdcrd' (protein+ligand + ten closest H2O + remaining bulk tip3p H2O) using a new topology created by ambpdb using one snapshot taken from out1.mdcrd. There is no problem viewing this trajectory containing all waters in box, the problem goes to strip water trajectory.
>
> I digress .... it is a fact that handling water molecules is challenging. I believe there is a solution to these issues.
>
> Amor
>
>
>
> --- On Fri, 20/8/10, Amor San Juan <amorsanjuan.yahoo.com> wrote:
>
> From: Amor San Juan <amorsanjuan.yahoo.com>
> Subject: PTRAJ: closest & strip in one script or separate
> To: "AMBER Mailing List" <amber.ambermd.org>
> Date: Friday, 20 August, 2010, 5:33 PM
>
>
> Is it possible to create a script in ptraj wherein both closest and strip commands are present?
>
> The script I used for closest:
>
> closest 10 :200
> image :WAT byres
> rms first mass :1-199
>
>
> The output trajectory (in order as it appears) contains :
> protein
> ligand
> waters(10) closest
> remaining tip3p water
>
> The script was able to find the ten closest waters near the ligand (:200) and write it in the output immediately afer the protein+ligand. The problem is the unselected waters (thousands and thousands) are still in the output. I tried to experiment on adding the strip command in above script but was not successful. Also, the other way around is from the output above I make a separate ptraj script to strip the waters
> 211-9199. Both are not successful.
>
> My aim is to get an output trajectory (P+L+W10), such that P=protein,L=ligand,W10=ten closest waters. And feed this into MMPBSA for binding energy calculation.
>
> Amor
>
>
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Aug 24 2010 - 05:00:05 PDT
