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