Re: PBC in lipid bilayer: receptor vs. ligand

From: Jeff Comer (jeffcomer_at_gmail.com)
Date: Thu Jan 28 2016 - 12:38:45 CST

You can have it all! I usually shift the protein to the origin and wrap
everything else around it. Try something like:

package require pbctools
set sel [atomselect top protein]
set all [atomselect top all]

set frameNum [molinfo top get numframes]
for {set f 0} {$f < $frameNum} {incr f} {
    molinfo top set frame $f

    # Shift the center to the origin.
    set cen [measure center $sel weight mass]
    $all moveby [vecinvert $cen]
}

# Wrap about the center.
pbc wrap -compound res -center origin -all
animate write dcd $outDcd beg 0 end -1 sel $all top

By translating your system and using various "pbc wrap" and "pbc unwrap"
commands (and sometimes "pbc join") you can do most anything.

Jeff

–––––––––––––––––––––––––––––––––––———————
Jeffrey Comer, PhD
Assistant Professor
Institute of Computational Comparative Medicine
Nanotechnology Innovation Center of Kansas State
Kansas State University
Office: P-213 Mosier Hall
Phone: 785-532-6311

On Thu, Jan 28, 2016 at 11:59 AM, R. Charbel MAROUN <
charbel.maroun_at_inserm.fr> wrote:

> Absolutely right. It wasn't a problem of writing or reading the dcd file.
> But, thing is, to be "fair" ;-), I'd have to do an unwrap for the lipids
> too. So:
>
>> set cell
>>
> {76.105133 86.859406 165.689102 90.000000 90.000000 90.000000}
>
>> pbc set {76.105133 86.859406 165.689102 90.000000 90.000000 90.000000}
>> -now
>> pbc unwrap -sel protein
>> pbc unwrap -sel "resname HME"
>> pbc unwrap -sel lipids
>>
>
> When doing that, protein, HME and lipid bi-layer get shifted in the plane
> of the membrane with respect to the water bi-layer. When unwrapping protein
> and ligand, both get shifted with respect to lipids and solvent. Guess
> can't have it all :-(.
>
> Cheers and thanx.
> ---
> R. Charbel MAROUN, Ph.D., H.D.R.
> UMR-S INSERM U1204/UEVE
> Structure et activité des biomolécules
> normales et pathologiques
> Université d'Evry-Val d'Essonne
> Bâtiment Maupertuis
> Rue du Père Jarlan
> 91000 EVRY
> FRANCE
> Tél: +33 1 69 47 76 64
> FAX: +33 1 69 47 02 19
> e-mail charbel.maroun_at_inserm.fr
>
> Le 28-01-2016 18:09, Josh Vermaas a écrit :
>
>> Exclude water from the wrapping (-sel "not water"). Unwrap works by
>> asking the following question for every atom: "Did you move more than
>> half a periodic box length from the previous frame?" If yes, it moves
>> the atom so that it didn't move more than half a box length in any
>> direction from the previous frame. Water moves very quickly, so
>> between two adjacent frames one of the atoms might have moved more
>> than half the box length, but the other two might not have. The
>> algorithm then dutifully moves the single atom, creating one or two
>> long bonds, which create the disk-like shape you are observing.
>>
>> -Josh
>>
>> On 01/28/2016 11:03 AM, R. Charbel MAROUN wrote:
>>
>>> Hi,
>>>
>>> Indeed, I played with the pbc parameters. I think I got the solution
>>> with "all" as the unwrap option:
>>>
>>> set cell
>>>>
>>> {76.105133 86.859406 165.689102 90.000000 90.000000 90.000000}
>>>
>>>> pbc set {76.105133 86.859406 165.689102 90.000000 90.000000 90.000000}
>>>> -now
>>>> pbc unwrap -sel all
>>>>
>>> 1.3% complete (frame 1)
>>> 100.0% complete (frame 149)
>>>
>>> That keeps receptor AND ligand together.
>>>
>>> After that I saved all atoms of all frames of the modified trajectory in
>>> dcd format using File Save Trajectory (Save all at once). But when reading
>>> the new trajectory back into VMD, I observe that the first frame is OK, but
>>> as time goes by, the water bi-layer literally "explodes" into an enormous
>>> disk-like bi-layer (I'm using the same psf file). Any clues? Another way to
>>> save a trajectory? (Maybe I should file this problem under another thread).
>>>
>>> Cheers.
>>>
>>> PS Yes, I had the whole trajectory loaded.
>>>
>>>
>>> ---
>>> R. Charbel MAROUN, Ph.D., H.D.R.
>>> UMR-S INSERM U1204/UEVE
>>> Structure et activité des biomolécules
>>> normales et pathologiques
>>> Université d'Evry-Val d'Essonne
>>> Bâtiment Maupertuis
>>> Rue du Père Jarlan
>>> 91000 EVRY
>>> FRANCE
>>> Tél: +33 1 69 47 76 64
>>> FAX: +33 1 69 47 02 19
>>> e-mail charbel.maroun_at_inserm.fr
>>>
>>> Le 28-01-2016 16:43, Josh Vermaas a écrit :
>>>
>>>> Just a quick question, do you have the whole trajectory loaded? I have
>>>> to admit, I'm a bit fuzzy on the algorithm of how unwrap works, but
>>>> there is a way of getting pbctools to do what you want. Play around
>>>> with it first (Norman's suggestion also looked good), since its much
>>>> cheaper to rewrap the trajectory than to run another simulation!
>>>> -Josh
>>>>
>>>> On 01/28/2016 05:26 AM, R. Charbel MAROUN wrote:
>>>>
>>>>> Hello namders,
>>>>>
>>>>> Following Josh's proposition, I used
>>>>> pbc unwrap -sel "resname HME"
>>>>> where HME is the name of the ligand.
>>>>>
>>>>> pbc is unwrapping HME from its present receptor-unpaired positions.
>>>>> That's leading in many cases to other receptor-unpaired positions rather
>>>>> than to the original receptor-paired positions. I guess this is due to the
>>>>> fact that pbc cannot know what was the original unwrapped position of the
>>>>> ligand.
>>>>>
>>>>> So, back to the first stage.
>>>>>
>>>>> Cheers,
>>>>> ---
>>>>> R. Charbel MAROUN, Ph.D., H.D.R.
>>>>> UMR-S INSERM U1204/UEVE
>>>>> Structure et activité des biomolécules
>>>>> normales et pathologiques
>>>>> Université d'Evry-Val d'Essonne
>>>>> Bâtiment Maupertuis
>>>>> Rue du Père Jarlan
>>>>> 91000 EVRY
>>>>> FRANCE
>>>>> Tél: +33 1 69 47 76 64
>>>>> FAX: +33 1 69 47 02 19
>>>>> e-mail charbel.maroun_at_inserm.fr
>>>>>
>>>>> Le 27-01-2016 16:22, Josh Vermaas a écrit :
>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> The simplest solution would be to turn wrapAll off, since what is
>>>>>> happening is that the center of mass of the segment (in this case just
>>>>>> the ligand) is going across the periodic boundary edge and NAMD is
>>>>>> wrapping the position around to the origin. However, since the
>>>>>> simulation has already been run, what I would recommend is to load it
>>>>>> in VMD, and unwrap just the ligand, and then save the trajectory
>>>>>> again. See the pbctools plugin
>>>>>> (http://www.ks.uiuc.edu/Research/vmd/plugins/pbctools/). This command
>>>>>> *should* bring the ligand back to where you expect it.
>>>>>>
>>>>>> pbc unwrap -sel "text to select the ligand in VMD"
>>>>>>
>>>>>> -Josh Vermaas
>>>>>>
>>>>>> On 01/27/2016 07:54 AM, R. Charbel MAROUN wrote:
>>>>>>
>>>>>>> Hello NAMD users:
>>>>>>>
>>>>>>> I have a receptor embedded in a lipid bilayer. The receptor has a
>>>>>>> ligand in its binding site. After several hundred ps of MD, the receptor
>>>>>>> (and the ligand) move away from the center of the membrane in the plane of
>>>>>>> the membrane. As I have applied PBC, when the receptor gets away far enough
>>>>>>> from the membrane, its image appears at the other side. Sometimes, it's the
>>>>>>> image of the ligand that appears in the other side. This depends on whether
>>>>>>> the receptor or the ligand attains the border of the cell. As a
>>>>>>> consequence, the receptor and ligand get unpaired for many frames, ie, the
>>>>>>> ligand is not in the binding site. Apparently, the ligand continues to
>>>>>>> behave as if in it. The problem comes if I want to measure, for ex.,
>>>>>>> distances between ligand and residues for those frames, or make statistics
>>>>>>> out of my trajectory. Below is the inp file I'm using.
>>>>>>>
>>>>>>> Is there any way around, such as imposing a constraint to the
>>>>>>> ligand, so that it follows the receptor, even if the ligand doesn't "hit"
>>>>>>> the border of the cell ?
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> structure C_3CAP-H2ODow-HME_sol_ion.psf
>>>>>>> coordinates C_3CAP-H2ODow-HME_sol_ion.pdb
>>>>>>>
>>>>>>> outputName step74_prod; # base name for output from this run
>>>>>>> # NAMD writes two files at
>>>>>>> the end, final coord and vel
>>>>>>> # in the format of
>>>>>>> first-dyn.coor and first-dyn.vel
>>>>>>>
>>>>>>> set inputname step73_prod;
>>>>>>> binCoordinates $inputname.restart.coor; # coordinates from
>>>>>>> last run (binary)
>>>>>>> binVelocities $inputname.restart.vel; # velocities from
>>>>>>> last run (binary)
>>>>>>> extendedSystem $inputname.restart.xsc; # cell dimensions
>>>>>>> from last run (binary)
>>>>>>>
>>>>>>> firsttimestep 2254880; # last step of previous run
>>>>>>> restartfreq 500; # 500 steps = every 1ps
>>>>>>> dcdfreq 1000;
>>>>>>> dcdUnitCell yes; # the file will contain unit
>>>>>>> cell info in the style of
>>>>>>> # charmm dcd files. if yes,
>>>>>>> the dcd files will contain
>>>>>>> # unit cell information in
>>>>>>> the style of charmm DCD files.
>>>>>>> xstFreq 500; # XSTFreq: control how often
>>>>>>> the extended systen configuration
>>>>>>> # will be appended to the XST
>>>>>>> file
>>>>>>> outputEnergies 125; # 125 steps = every 0.25ps
>>>>>>> # The number of timesteps
>>>>>>> between each energy output of NAMD
>>>>>>> outputTiming 500; # The number of timesteps
>>>>>>> between each timing output shows
>>>>>>> # time per step and time to
>>>>>>> completion
>>>>>>>
>>>>>>> # Force-Field Parameters
>>>>>>> paraTypeCharmm on; # We're using charmm type
>>>>>>> parameter file(s)
>>>>>>> # multiple definitions may be
>>>>>>> used but only one file per definition
>>>>>>>
>>>>>>> # parameters
>>>>>>> /usr/local/vmd-.9/lib/vmd/plugins/noarch/tcl/readcharmmpar1.2/par_all27_prot_lipid_na.inp
>>>>>>>
>>>>>>> parameters
>>>>>>> /home/cmaroun/toppar/toppar_27/par_all27_prot_lipid_cholesterol_HME_TIP3.prm
>>>>>>> parameters /home/cmaroun/readcharmmpar1.2/par_all36_lipid.prm
>>>>>>>
>>>>>>> # These are specified by CHARMM
>>>>>>> exclude scaled1-4 # non-bonded exclusion policy
>>>>>>> to use "none,1-2,1-3,1-4,or scaled1-4"
>>>>>>> # 1-2: all atoms pairs that
>>>>>>> are bonded are going to be ignored
>>>>>>> # 1-3: 3 consecutively bonded
>>>>>>> are excluded
>>>>>>> # scaled1-4: include all the
>>>>>>> 1-3, and modified 1-4 interactions
>>>>>>> # electrostatic scaled by
>>>>>>> 1-4scaling factor 1.0
>>>>>>> # vdW special 1-4 parameters
>>>>>>> in charmm parameter file.
>>>>>>> 1-4scaling 1.0
>>>>>>> switching on
>>>>>>> vdwForceSwitching yes; # New option for force-based
>>>>>>> switching of vdW
>>>>>>> # if both switching and
>>>>>>> vdwForceSwitching are on CHARMM force
>>>>>>> # switching is used for vdW
>>>>>>> forces.
>>>>>>>
>>>>>>> # You have some freedom choosing the cutoff
>>>>>>> cutoff 12.0; # may use smaller, maybe 10.,
>>>>>>> with PME
>>>>>>> switchdist 10.0; # cutoff - 2.
>>>>>>> # switchdist - where you
>>>>>>> start to switch
>>>>>>> # cutoff - where you stop
>>>>>>> accounting for nonbond interactions.
>>>>>>> # correspondence in charmm:
>>>>>>> # (cutnb,ctofnb,ctonnb =
>>>>>>> pairlistdist,cutoff,switchdist)
>>>>>>> pairlistdist 14.0; # stores the all the pairs
>>>>>>> with in the distance it should be larger
>>>>>>> # than cutoff( + 2.)
>>>>>>> stepspercycle 20; # 20 redo pairlists every ten
>>>>>>> steps
>>>>>>> pairlistsPerCycle 2; # 2 is the default
>>>>>>> # cycle represents the number
>>>>>>> of steps between atom reassignments
>>>>>>> # this means every 20/2=10
>>>>>>> steps the pairlist will be updated
>>>>>>>
>>>>>>> # Integrator Parameters
>>>>>>> timestep 2.0; # fs/step
>>>>>>> rigidBonds all; # Bound constraint all bonds
>>>>>>> involving H are fixed in length
>>>>>>> nonbondedFreq 1; # nonbonded forces every step
>>>>>>> fullElectFrequency 1; # PME every step
>>>>>>>
>>>>>>> wrapWater on; # wrap water to central cell
>>>>>>> wrapAll on; # wrap other molecules too
>>>>>>> wrapNearest off; # use for non-rectangular
>>>>>>> cells (wrap to the nearest image)
>>>>>>>
>>>>>>> # PME (for full-system periodic electrostatics)
>>>>>>> source checkfft.str
>>>>>>> margin 2.5;
>>>>>>> PME yes;
>>>>>>> PMEInterpOrder 6; # interpolation order (spline
>>>>>>> order 6 in charmm)
>>>>>>> PMEGridSizeX $fftx; # should be close to the cell
>>>>>>> size
>>>>>>> PMEGridSizeY $ffty; # corresponds to the charmm
>>>>>>> input fftx/y/z
>>>>>>> PMEGridSizeZ $fftz;
>>>>>>>
>>>>>>> # Constant Temperature Control
>>>>>>> set temp 323.15;
>>>>>>> langevin on; # langevin dynamics
>>>>>>> langevinDamping 1.0; # damping coefficient of 1/ps
>>>>>>> (keep low)
>>>>>>> langevinTemp $temp; # random noise at this level
>>>>>>> langevinHydrogen no; # don't couple bath to
>>>>>>> hydrogens
>>>>>>> reinitvels $temp;
>>>>>>>
>>>>>>> # Constant Pressure Control (variable volume)
>>>>>>> useGroupPressure yes; # use a hydrogen-group based
>>>>>>> pseudo-molecular viral to calcualte pressure and
>>>>>>> # has less fluctuation, is
>>>>>>> needed for rigid bonds (rigidBonds/SHAKE)
>>>>>>> useFlexibleCell yes; # yes for anisotropic system
>>>>>>> like membrane
>>>>>>> useConstantRatio yes; # keeps the ratio of the unit
>>>>>>> cell in the x-y plane constant A=B
>>>>>>>
>>>>>>> langevinPiston on; # Nose-Hoover Langevin piston
>>>>>>> pressure control
>>>>>>> langevinPistonTarget 1.0; # target pressure in bar 1atm
>>>>>>> = 1.01325bar
>>>>>>> langevinPistonPeriod 50.0; # oscillation period in fs.
>>>>>>> correspond to pgamma T=50fs=0.05ps
>>>>>>> # f=1/T=20.0(pgamma)
>>>>>>> langevinPistonDecay 25.0; # oscillation decay time.
>>>>>>> smaller value correspons to larger random
>>>>>>> # forces and increased
>>>>>>> coupling to the Langevin temp bath.
>>>>>>> # Equall or smaller than
>>>>>>> piston period
>>>>>>> langevinPistonTemp $temp; # coupled to heat bath
>>>>>>> run 1500000; # 3ns
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>

This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:21:47 CST