[AMBER] Creating new topology file from close.mdcrd file

From: Amor San Juan <amorsanjuan.yahoo.com>
Date: Tue, 17 Aug 2010 02:42:40 -0700 (PDT)

Thanks Jason for always giving time and attention.

The input script I used for closest is:

----
trajin CLOin.mdcrd
trajout CLOout.mdcrd
center :1-200 mass origin
image origin center
solvent byres :WAT
closest 10 :200 first
-----
The ligand numbering is 200, where I want 10 closest water around it. I got the output trajectory but I bumped into another problem of generating the new topology file corresponding to "CLOout.mdcrd". This new topology file shall be obtained from a new pdb file (P+L+W10) where W10 is 10 water molecules. I wanted to do the ff steps:
1.Take one snapshot from "CLOout.mdcrd" by extraction of trajectory using mmpbsa.pl. Convert snapshot to pdb using ambpdb.
2. Modify the pdb file such that only 10 water molecules remained, together with prot+lig
3.Create new topology using the pdb generated in step2.
File CLOin.mdcrd input trajectory has 37581879 bytes while CLOout.mdcrd output trajectory has only 4161681 bytes. Input traj contains 9000 TIP3P model waters, and after executing closest command the output is expected to contain only 10 closest water molecules (but Im not sure of this assumption). 
My question is: does the output trajectory from "closest" command gives only the specified number of water? Or, the selected 10 closest water are written immediately after prot+lig, and followed by the rest of thousands of tip3p models which are farthest. I need to know precisely so that I will track down the error of generating snaphots from CLOout.mdcrd.
In this case, the CLOin.mdcrd has (Prot+Lig+Wat)=PLW.prmtop contains 30324 atoms. This PLW.prmtop does not correspond to CLOout.mdcrd, because whenever I tried generating snapshots there is mismatch of atoms. I coudnt be able to generate a pdb file from step1 above.
Mismatch atoms in generating pdb seems to me that PLW10.prmtop should be used, but to where I could get it as this is the reason I am doing steps 1-3 to get it.
Please anybody help me get from this bump. I have exhausted the entire day.
Amor
--- On Mon, 16/8/10, Jason Swails <jason.swails.gmail.com> wrote:
From: Jason Swails <jason.swails.gmail.com>
Subject: Re: [AMBER] How to identify water molecules at the interface ? follow up ....
To: "AMBER Mailing List" <amber.ambermd.org>
Date: Monday, 16 August, 2010, 9:03 PM
Hello,
On Mon, Aug 16, 2010 at 2:19 AM, Amor San Juan <amorsanjuan.yahoo.com>wrote:
> I want to write a follow-up in regard to my previous post below.
>
> Aside from ptraj hbond analysis, I also tried several ptraj commands
> including:
> a) grid = the aim of this command is to generate water density at hydration
> b)watershell= the aim of this command is to count the number of water
> molecules at 1st/2nd solvation shell
> c) closest= the aim of this command is to quantify the number of water
> molecules at 3.5A and 5A.
>
> In all those commands above, I am unable to identify the numbering of water
> buried at interface. Although the closest command can quantify the water
> molecules at specific distance, its numbering of water molecules were
> altered compared to the original numbering.
>
> My aim is to get the identity of the water molecules buried between the
> ligand-protein. Once the water identity is found, I will use this as input
> together with prot+lig for binding energy calculation in mmpbsa. What I mean
> to say in water identity is simply the water identity as reflected on its
> atom/residue number.
>
You don't need this really.  Just create a new prmtop and trajectory file
with only the proper number of waters that you kept with the "closest"
command.  Then it doesn't matter what the original numbering was in the
prmtop/trajectory file(s), it's simply the new numbers it has in the new
trajectory.  Then feed this into MMPBSA and don't tell it to strip any
waters/ions.  For example, if you're using MMPBSA.py, set strip_mdcrd=0 to
avoid stripping the trajectory file of any waters/ions.  For mm_pbsa.pl,
just include the water atoms as part of the ligand or receptor atoms.
Hope this helps,
Jason
> amor
>
>
> --- On Mon, 16/8/10, Amor San Juan <amorsanjuan.yahoo.com> wrote:
>
> From: Amor San Juan <amorsanjuan.yahoo.com>
> Subject: How to identify water molecules at the interface ?
> To: "AMBER Mailing List" <amber.ambermd.org>
> Date: Monday, 16 August, 2010, 11:51 AM
>
> I have read several posts regarding water calculations in several ways for
> different purposes done. But it seems to me that the goal I want to achieve
> needs relevant feedback from the forum.
>
> I have a system (protein + peptide + ligand) with 200 residues in total and
> note that ligand has a phosphate group. I want to identify the water
> molecules (TIP3P models) found at the interface between the ligand and
> protein. I need to identify which water molecules are at the interface so I
> can use it as input for binding energy calculation by MMPBSA. Here are the
> steps formulated for water identification:
>
> 1) Calculate the water within 3.5Angstroms using the trajectory input. The
> expected outcome is (P+L+Wi) such that P=protein, L=ligand, and Wi=water at
> ith frame. From the output, remove all hydrogen molecules in water.
>
> 2)Calculate accessibility by
>  filtering out WAT.O (Water oxygen) within >0.50Angstroms
>
> 3) Repeat step 2 to get only the buried water molecules remaining at
> interface.
>
> One suggestion I got is to make a script to do in one shot the above three
> steps. I think about it but am not sure if one straightforward script can do
> the task.
>
> I approached the problem by trying to do it one-step-at-a-time. For step 1
> above, I implemented ptraj. Using ptraj h-bond distance 3.5 analysis, it
> gave me an output of P+L+Wi, however looking at Wi seems to be so enormous
> from the list as expected.
> For step2, I need the trajectory that contains P+L+Wi with stripped
> hydrogens generated from step1. I have not proceeded to step2 yet.
>
> To get the trajectory P+L+Wi, it seems the way out is to use the
> ptraj-hbond output to select interacting atoms, but as noted Wi in the list
> is quite long to do manual selection one-by-one.
>
>  I also thought of using LIGPLOT program to
>  identify buried water molecules but my input is a trajectory not a pdb. I
> know I can generate pdb from snapshots in trajectory but how would you take
> it for several thousands of snapshots pdb? Bottomline, I only need to
> identify the water at interface based on the trajectory.
>
> Amor
>
>
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
-- 
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Graduate Student
352-392-4032
_______________________________________________
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 17 2010 - 03:00:03 PDT
Custom Search