PBC in lipid bilayer: receptor vs. ligand

From: R. Charbel MAROUN (charbel.maroun_at_inserm.fr)
Date: Wed Jan 27 2016 - 07:54:34 CST

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

-- 
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

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:20:09 CST