WRITTING MULTI_FRAME PDB FILE FOR DYNAMIC SELECTION TRAJECTORY

133 views
Skip to first unread message

Mireia Segado

unread,
Oct 4, 2023, 4:36:10 AM10/4/23
to MDnalysis discussion
Dears,

I would like to select the first solvation shell to do a specific  analysis generating a multiframe -pdb.. But the pdb generated the waters selected are relate to first frame. Which is the best manner to do it?

B12F12 = mda.Universe(file.pdb)
B12F12_cluster = B12F12.select_atoms("resname BF1",updating= True)
B12F12_1stSS = B12F12.select_atoms("resname SOL and around 6.0 resname BF1", updating = True)

with mda.Writer("B12F12_1stShell.pdb", multiframe = True) as W:
    for ts in B12F12.trajectory:
        W.write(B12F12_1stSS)

Thanks in advance

Mire

Mireia Segado

unread,
Oct 4, 2023, 8:43:20 AM10/4/23
to MDnalysis discussion
Dears, 

I also tried reading .xtc and .tpr.

I would like to clarify: The problem of writing selected coordinates to the pdb seems correct, the problem is that in each frame-MODEL in the generated PDB there are not the Dimensions of the simulation box for each MODEL. The dimensions of the box are only written for model1. Its a NPT trajectory.

Is there any solution to solve this problem and write the box dimensions for each frame?

Thank you!

Mire

Oliver Beckstein

unread,
Oct 5, 2023, 7:56:44 AM10/5/23
to mdnalysis-discussion
Hello Mire,

MDAnalysis cannot deal with systems with varying numbers of atoms easily. In particular, we can only write MODEL entries with the same number of atoms. I am not 100% sure if the PDB standard allows MODEL entries with different numbers of atoms but if that’s what you want  I can think of two ways to somehow get close to what you want. Both appraoches require extra work:

1. Write a single PDB file for each frame and then manually combine the files into a multiframe PDB.
2. Write a multiframe PDB file with ALL atoms in the system but use the occupancy field to indicate which atoms are selected.

For 1 (which will be quite slow):

B12F12 = mda.Universe(“file.pdb")
B12F12_1stSS = B12F12.select_atoms("resname SOL and around 6.0 resname BF1", updating=True)

for ts in B12F12.trajectory:
    filename = f"B12F12_1stShell_{ts.frame:06d}.pdb”
    B12F12_1stSS.write(filename)    # can directly use AtomGroup.write()

Now you have to write your own code to concatenate the PDB files and insert the appropriate MODEL/ENDMDL records. Eg in the shell something like the following (NOT TESTED!!)

grep “^CRYST1” B12F12_1stShell_000001.pdb > multi.pdb
for PDB in *.pdb; do
(echo “MODEL”; grep “^ATOM” $PDB; echo “ENDMDL”) >> multi.pdb
        done

or probably better to do this in Python. In any case, you have to figure out how to combine the individual PDB files.

For 2 (this *should* work although I wasn’t able to try it out):

B12F12 = mda.Universe(“file.pdb”)
B12F12.add_TopologyAttr(“occupancies”)


B12F12_1stSS = B12F12.select_atoms("resname SOL and around 6.0 resname BF1", updating=True)

with mda.Writer("B12F12_1stShell.pdb", multiframe=True) as W:
   for ts in B12F12.trajectory:
       # set everything to “excluded”
       B12F12.atoms.occupancies = 0
       
       # set solvation shell to included
       B12F12_1stSS.occupancies = 1

       # write whole system frame with occupancies
       W.write(B12F12.atoms)

Let us know if this solves your problem.

Oliver

-- 
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/9cb3bb03-5a3f-44f3-8bb0-7f806210ca02n%40googlegroups.com.
--
Oliver Beckstein (he/his/him)

email: orbe...@mdanalysis.org
twitter: @orbeckst
GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project
Reply all
Reply to author
Forward
0 new messages