From: R. Charbel MAROUN (charbel.maroun_at_inserm.fr)
Date: Thu Jan 28 2016 - 11:59:14 CST
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