CHARMM Development Project
Posted By: blubbi Selecting only lipids under protein - 09/14/10 03:07 PM
Hi Charmm users,

I tried to come up with a smart selection to only select the lipids under my protein. Since the protein moves on the X/Y-Plane I had the following plan:

Get the center of mass for my protein, extract the X/Y-Coordinates and select all lipids within 15A of these Coordinates. Is this possible within one selection?

I am interested in the dynamics, density, overall charge, etc.. of the lipids below the protein compared to the lipids outside the selection.

The selection would have to be updated every step. Is this possible with a selection or do I have to iterate every step and save which lipids are at which step under my protein and then calculate ${WHATEVER} separately for this specific residue at the specific steps?

I hope I could make myself clear.

Thanks for any suggestions.
Cheers
Bjoern
Posted By: lennart Re: Selecting only lipids under protein - 09/14/10 03:29 PM
Something like this snippet should do what you are asking about (but I am not sure that your criterion for lipids being under the protein is optimal):

label loop
traj read
coor stat sele segid prot end
define lipunder sele segi lipi .and. point ?XAVE ?YAVE ?ZAVE CUT 15 end
.
.
.
goto loop
Posted By: rmv Re: Selecting only lipids under protein - 09/14/10 05:07 PM
For time series analysis involving the notion of nearby molecules, the shell feature should be useful; see shell.doc and correl.doc for more info.
Posted By: lennart Re: Selecting only lipids under protein - 09/14/10 06:13 PM
I agree, for many analyses of solvent behavior around a solute the existing tools (COOR ANALysis or SHELL) do a very good job. For the kinds of analyses indicated by the original poster a CHARMM loop is likely to be the better solution; in particular since the concept of "under" is not so clearly tied to "being close to".
Posted By: rmv Re: Selecting only lipids under protein - 09/14/10 08:52 PM
True, but in many contexts, such as the above snippet, atom selection is not always "image aware" (shell code is), which may be an issue. Depending on the cell size, a 15 A shell based on lipid P and N atoms (assuming PC lipid) may be fairly reasonable.

For a loop over frames, the image cell boundary problem could also be handled by a coordinate translation in the xy plane followed by an image update, but that would slow things down a bit more.

I often try several different approaches in this situation (something to do with eggs and baskets smile ).
Posted By: blubbi Re: Selecting only lipids under protein - 09/15/10 09:05 AM
Originally Posted By: lennart

label loop
traj read
coor stat sele segid prot end
define lipunder sele segi lipi .and. point ?XAVE ?YAVE ?ZAVE CUT 15 end
goto loop


Looks promising, thanks a lot!
I would have to add an offset to ZAVE or increase the CUT since the max diameter in the current selection will be somewhat over the membrane interface.


Originally Posted By: lennart
I agree, for many analyses of solvent behavior around a solute the existing tools (COOR ANALysis or SHELL) do a very good job. For the kinds of analyses indicated by the original poster a CHARMM loop is likely to be the better solution; in particular since the concept of "under" is not so clearly tied to "being close to".


This is sadly so true, since the protein has a "hole" in the middle an if I would use nearby those lipids under the center of protein would most likely be excluded. Maybe this ASCII-Art can clear things up:

Code:
       _______
      /  ___  \  Protein
     /__/   \__\
||||||||||||||||||| Membrane
|||||||||||||||||||
      |-------|
   Lipids of interest


Originally Posted By: rmv
True, but in many contexts, such as the above snippet, atom selection is not always "image aware" (shell code is), which may be an issue. Depending on the cell size, a 15 A shell based on lipid P and N atoms (assuming PC lipid) may be fairly reasonable.

For a loop over frames, the image cell boundary problem could also be handled by a coordinate translation in the xy plane followed by an image update, but that would slow things down a bit more.

I often try several different approaches in this situation (something to do with eggs and baskets smile ).


Thanks for pointing out this problem.
Could you give an code example for"coordinate translation in the xy plane". I know the purpose behind the translation and in fact I probably would have run into this pitfall, since my protein is is moving around and reaches the end of the box here and there so I guess I have to deal with the boundary conditions. crazy

How much would this slow things down? It's going to be a very long simulation (so far ~400ns but still crunching to reach some microseconds and I am dreaming of milliseconds)
The simulations are consecutive runs of 1ns slices. So for each unreduced trajectory (1ns, 500000 steps, write out freq. 500, integration time 2 femtoseconds) has the size of 1.5G.

Sadly, running multiple small simulations is not an option here cry

If the slowdown is way to big I maybe should go with a more inaccurate selection (COOR ANALysis or SHELL) in favor for a reasonable time.

I think I'll play a bit around with both options and see what happens.

Thanks a lot!
Bjoern
Posted By: lennart Re: Selecting only lipids under protein - 09/15/10 02:12 PM
You ought to be able to use MERGE RECEnter to get the protein in the middle of the box in all frames.
To save diskspace you could increase the write out frequency, say to 2500. This still gives you 200 frames/ns which is plenty for such a long simulation.
Posted By: blubbi Re: Selecting only lipids under protein - 09/15/10 03:05 PM
Originally Posted By: lennart
You ought to be able to use MERGE RECEnter to get the protein in the middle of the box in all frames.

I am not entirely familiar how recenter in detail works, but I only used it to eliminates translation of the entire system. But since the Protein is moving relative to the membrane this would have no effect on the protein moving around on the membrane. Please correct me if understood recentering wrong. Another concern I do have regarding the recentering on the protein. If I should choose to recenter on the protein will this not influence the lipids I see passing/resting/following below the protein, will it? I am confused and I couldn't dissipate my doubt by reading the RECEnter notes in dynam.doc.

Originally Posted By: lennart

To save diskspace you could increase the write out frequency, say to 2500. This still gives you 200 frames/ns which is plenty for such a long simulation.

Well, I am not sure about this, since I once ran into trouble with a too low resolution (1ns) while calculating relaxation times. Comparing the simulation to the same setup with 100ps time resolution yielded different results. So far we didn't investigate any further but the difference might be related to the time resolution. So I'd rather throw in some more discs then throw away time by recalculating with a higher time resolution.

Deleting trajectories and storing reduced ones is quick and always possible ;-)

Thanks for your help and suggestions.
I appreciate it very much!

Cheers
Bjoern
Posted By: rmv Re: Selecting only lipids under protein - 09/15/10 05:06 PM
I hesitate to use MERGE COOR RECENTER, since the the protein will be shifted in z as well as x and y, which may not be desirable. Also, if (1) the bilayer remains centered around z=0 during the simulation, and (2) the distance between the protein an bilayer surface is variable, you may wish to use a fixed point in z instead of ?ZAVE, which is the center of the protein. If the protein is always above (+z) the bilayer, ?ZMIN is a better choice as well. For the x,y shift--

traj read
coor stat sele segid prot end
coor tran xdir -?XAVE ydir -?YAVE
update imgfrq 1 ...
define lipunder sele segi lipi .and. point 0. 0. ?ZAVE CUT 15 end
Posted By: blubbi Re: Selecting only lipids under protein - 09/16/10 01:25 PM
Well, this looks like a plan.

Thanks a lot for your suggestions!

Cheers
Bjoern
Posted By: rmv Re: Selecting only lipids under protein - 09/16/10 04:43 PM
If you intend to do a number of different analyses using this approach, it might be a good idea to create new trajectories with the x,y centering already performed. That way you can eliminate the overhead of the image update from subsequent runs.
Posted By: blubbi Re: Selecting only lipids under protein - 09/17/10 12:17 PM
Originally Posted By: rmv
If you intend to do a number of different analyses using this approach, it might be a good idea to create new trajectories with the x,y centering already performed. That way you can eliminate the overhead of the image update from subsequent runs.


Good point! Thanks a lot!

Cheers
Bjoern
Posted By: blubbi Re: Selecting only lipids under protein - 01/19/11 10:47 AM
Hi again.

I have a question regarding the "update" command.
As long as I am not interested in nonbonded properties like H-Bonds and other electrostatic stuff is it safe to ignore the update command?

Lets say I read the trajectory without setting up the crystal information, iterate it step by step while re-centering every step and then selecting the Lipids under my protein.

Is it necessary to use the "update" command?
And is my assumption right that the update command only updates the nonbonded properties (Which would make no sense when there is no crystal information)

Thanks a lot,
Bjoern
Posted By: lennart Re: Selecting only lipids under protein - 01/19/11 10:54 AM
If you are only doing geometry related analyses I don't see any need for the update command. COOR HBONd also works fine without UPDAte.
Posted By: blubbi Re: Selecting only lipids under protein - 01/19/11 11:40 AM
Thanks for confirming my assumption.

Cheers,
Bjoern
Posted By: rmv Re: Selecting only lipids under protein - 01/19/11 03:01 PM
That is not correct if you intend to look at protein:lipid interactions, unless you have a new trajectory created with MERGE COOR RECENTER. Regardless of whether the analysis is geometric or energetic, if you move the system to center the protein, if it was protruding from the unit cell, there will be many missing lipids UNLESS you do an image update to repack the lipids. How do you intend to recenter w/o the crystal information?

Also, COOR HBOND is only image aware if if you've set up images.
Posted By: blubbi Re: Selecting only lipids under protein - 01/19/11 05:18 PM
Thanks a lot for the rectification.

One thing I am curious about:
When running a simulation for a long time it might grow or shrink over time so the crystal data (side length of box in particular) I used to setup the simulation is "outdated".

When using "correl" or "update" with crystal information I ran into the following error:
Code:
      ***** LEVEL -1 WARNING FROM <TSTCRD> *****
      ***** SOME ATOM COORDINATES UNDEFINED OR OUT OF RANGE
      ******************************************
      BOMLEV (  0) IS REACHED - TERMINATING. WRNLEV IS  5

Here's the example code:
Code:
! Read in Topology and  Parameter files
stream ../read_psf_and_toppar.str

! Setup Crystal
crystal define TETR 93.7831 93.7831 144.964 90 90 90
open unit 1 read card name crystal.xtl
crystal read unit 1 card
close unit 1

! Image centering
image byresid XCEN 0.0 YCEN 0.0 ZCEN 0.0 sele ALL .and. .not. segid PRO* .and. .not. segid RNAA end
image byseg XCEN 0.0 YCEN 0.0 ZCEN 0.0 sele segid RNAA .or. segid PRO* end

open unit 99 file read name concatenated_recentered.dcd
traj query unit 99

calc NSTEP ?nstep
set START ?start
calc STOP ?nstep
set SKIP ?skip

calc mxa 9164 + 40334
calc mxs 6
calc mxt @NSTEP

correl maxtim @MXT maxa @MXA maxs @MXS

define MS sele segid MEMB end
define PS sele segid PRO* end

end
stop


The crystal.xtl was written with the "write crystal" command after solvating, neutralizing, centering and shaping the system.

When I omit the "read crystal.xtl" line, the error is gone.

Adding the "noupdate" option to "correl" omits the error as expected.

Or does this not at all relate to the side length of the box?

Cheers,
Bjoern
Posted By: lennart Re: Selecting only lipids under protein - 01/19/11 05:30 PM
This is completely unrelated to your box and your trajectory. correl can be used for energy related analyses, and by default calls the update routine to have various lists filled. Since in your case there are no coordinates present the error message you see is issued. You correctly added the noupdate keyword to avoid this; you could also have read in coordinates prior to entering the correl module, but the noupdate method is faster and cleaner.
Posted By: blubbi Re: Selecting only lipids under protein - 01/19/11 05:59 PM
While checking my other scripts I realized my mistake.

I simply missed to read in the initial coordinates as I do usually.

Your answer came to fast.
Sorry for the dumb mistake and bothering the forum with it.

But it's good to know that adding the "noupdate" option is the more appropriate way.

Thanks a lot - again -
Bjoern
© CHARMM forums