Re: PBC in lipid bilayer: receptor vs. ligand

From: R. Charbel MAROUN (charbel.maroun_at_inserm.fr)
Date: Thu Jan 28 2016 - 05:26:33 CST

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