Thank you Jason for your reply.
I'm running some test simulations right now, to find out, where is the
problem.
First of all - I've changed the prepration script from tleap, so I'm not
using the combine command. I'm thinking that the combine command is sure
messing something with my system. When I prepare the system without using
the combine command, density and volume are equilibrating quicker.
It looks like I don't quite understand, how does the iwrap work. I've
looked into my .crd files and I've found out, that the water box isn't
getting bigger at all.
$ tail -1 protein_complex.inpcrd
74.2459799 74.2459799 74.2459799 109.4712190 109.4712190 109.4712190
$ tail -1 emin3.rst (minimization)
74.2459799 74.2459799 74.2459799 109.4712190 109.4712190 109.4712190
$ tail -1 emd1.crd (heating)
74.2459799 74.2459799 74.2459799 109.4712190 109.4712190 109.4712190
$ tail -1 emd2.crd (equlibration)
71.2807594 71.2807594 71.2807594 109.4712190 109.4712190 109.4712190
$ tail -1 prod1.crd (production)
71.3345282 71.3345282 71.3345282 109.4712190 109.4712190 109.4712190
These values are the dimensions of my water box, am I right?
So it isn't getting any bigger during simulation.
But, as I can see in my prod1.out file, there is something like that:
...
wrapping first mol.: -713.94378 -0.00003 -466.34613
wrapping first mol.: -714.26667 -0.00003 -466.55705
wrapping first mol.: -714.26667 -0.00003 -466.55705
wrapping first mol.: -713.69154 -0.00003 -466.18137
wrapping first mol.: -713.69154 -0.00003 -466.18137
wrapping first mol.: -713.29150 -0.00003 -465.92006
wrapping first mol.: -713.29150 -0.00003 -465.92006
wrapping first mol.: -713.46314 -0.00003 -466.03218
wrapping first mol.: -713.46314 -0.00003 -466.03218
wrapping first mol.: -713.77968 -0.00003 -466.23894
wrapping first mol.: -713.77968 -0.00003 -466.23894
wrapping first mol.: -713.47238 -0.00003 -466.03822
wrapping first mol.: -713.47238 -0.00003 -466.03822
wrapping first mol.: -713.80941 -0.00003 -466.25836
wrapping first mol.: -713.80941 -0.00003 -466.25836
wrapping first mol.: -713.92023 -0.00003 -466.33075
wrapping first mol.: -713.92023 -0.00003 -466.33075
wrapping first mol.: -690.25657 33.66100 -524.72372
wrapping first mol.: -690.25657 33.66100 -524.72372
wrapping first mol.: -690.16879 33.65672 -524.65698
wrapping first mol.: -690.16879 33.65672 -524.65698
wrapping first mol.: -689.55424 33.62675 -524.18981
wrapping first mol.: -689.55424 33.62675 -524.18981
wrapping first mol.: -689.39230 33.61886 -524.06671
wrapping first mol.: -689.39230 33.61886 -524.06671
wrapping first mol.: -689.54014 33.62607 -524.17910
wrapping first mol.: -689.54014 33.62607 -524.17910
wrapping first mol.: -689.78757 33.63813 -524.36719
wrapping first mol.: -689.78757 33.63813 -524.36719
wrapping first mol.: -689.75031 33.63632 -524.33886
wrapping first mol.: -689.75031 33.63632 -524.33886
wrapping first mol.: -689.30142 33.61442 -523.99762
wrapping first mol.: -689.30142 33.61442 -523.99762
wrapping first mol.: -689.83884 33.64063 -524.40616
wrapping first mol.: -689.83884 33.64063 -524.40616
wrapping first mol.: -690.40235 33.66811 -524.83453
wrapping first mol.: -690.40235 33.66811 -524.83453
wrapping first mol.: -690.27845 33.66207 -524.74034
wrapping first mol.: -690.27845 33.66207 -524.74034
wrapping first mol.: -690.22461 33.65945 -524.69942
wrapping first mol.: -690.22461 33.65945 -524.69942
wrapping first mol.: -690.39049 33.66753 -524.82552
wrapping first mol.: -690.39049 33.66753 -524.82552
wrapping first mol.: -690.18416 33.65747 -524.66867
wrapping first mol.: -690.18416 33.65747 -524.66867
wrapping first mol.: -690.37099 33.66658 -524.81069
...
Can you tell me, what is going on here?
Is it something related with the periodic boundaries?
The volume is OK, the density is now OK, only the potential and total
energy are decreasing.
I thought that something is going on with the water box, because the
density was growing, volume was decreasing, and total energy was decreasing
too. And after 50ns the system was barely moving. I thought before - that
this is because there was no solvent so there were no interactions.
Right now I don't have a clue, what is going on or how to fix it.
Best regards,
Karolina
2015-10-21 14:10 GMT+02:00 Jason Swails <jason.swails.gmail.com>:
> On Wed, Oct 21, 2015 at 5:13 AM, Karolina Markowska <
> markowska.kar.gmail.com
> > wrote:
>
> > 2015-10-19 15:16 GMT+02:00 David A Case <dacase.scarletmail.rutgers.edu
> >:
> >
> > > On Mon, Oct 19, 2015, Karolina Markowska wrote:
> > > >
> > > > I have a problem with my system. It's a small globular protein (286
> > aa).
> > > > When I was doing short simulations (10 ns) there was no problem, but
> > > right
> > > > now, when I want to do 50 ns simulations (5 x 10 ns) I see there is
> > > > something wrong with the system.
> > > > The volume is getting smaller, density is increasing and the total
> > energy
> > > > is decreasing very fast.
> > >
> > > How much (numerically) is the density increasing? Are you starting the
> > > 50ns
> > > simulations from the end of the 10ns run? Are you using sander or
> pmemd?
> > > What do you see when you visualize the trajectory? Look *really*
> > > carefully to
> > > see if there are any differences between the 10ns simulations and the
> > > 5x10ns
> > > simulations.
> > >
> > > Now I see that there is also something funny with the density. It keep
> > increasing from 1.03 up to 1.053 and then it decrease.
> > I'm starting a new run. I'm not using the old 10ns simulation files. I'm
> > using pmemd.CUDA.
> >
>
> Why is this surprising? It's not unusual for system densities to be this
> high for protein systems (or even higher).
>
>
> Well, when I use this script:
> >
> > reference protein_complex.inpcrd
> > trajin prod1.traj.gz
> > trajin prod2.traj.gz
> > trajin prod3.traj.gz
> > trajin prod4.traj.gz
> > trajin prod5.traj.gz
> > autoimage
> > center :1-286 mass origin
> > image origin center familiar
> >
>
> These last two commands, "center" and "image" are probably redundant.
> "autoimage" usually does what you want it to.
>
>
>
> > rms ToRef :1-286.CA,C,N= reference out rmsd.arg mass
> > atomicfluct out rmsf.arg :1-286 byres
> > trajout protein_f.binpos binpos
> > go
> >
> > everything looks fine.
> > But when I get rid of autoimage, center and image commands, I see that my
> > water box is spinning over the center of the protein and the protein is
> > sticking out of the water box.
> >
>
> This is an imaging artifact, as proven by the fact that "autoimage"
> removes this problem. It has no effect on the simulation itself.
>
>
>
> > I think that this spinning octahedral could be something funny.
> >
>
> I really don't know what this means... it doesn't make sense to me that
> autoimage could remove box rotations (the unit cell should not rotate --
> the unit cell vectors should be fixed in direction through the whole
> simulation, and only the lengths allowed to change if using a barostat).
> Unless I'm missing something...
>
>
>
> > But probably this is not the main problem. The problem is, I think, that
> I
> > probably have free spaces between water from 3D-RISM and tleap waterbox.
> > I can see them, whet I make a pdb file from my .prmtop and .crd files
> from
> > each step of simulation. After minimization I've got this free spaces,
> but
> > during heating and equilibration everything looks OK to me.
>
>
> The barostat will press out these empty spaces, and it should be fine. In
> fact, tleap puts a lot of space between the edge of the box and the closest
> atoms, so there is going to be fairly large gaps if you visualize the
> system replicated in all directions. This leads to a *very* low internal
> pressure and simulation boiling if you run at constant volume. But using a
> Barostat will decrease the volume of the system until you force out these
> vacuum bubbles.
>
>
>
> > And when I try
> > to make a .pdb file from production coordinates file, I see that the
> > waterbox looks quite normal, except it have a hole on one side, and my
> > protein sticking out form another side of the box.
> > That's kinda weird.
> >
>
> Again, it's an imaging artifact.
>
>
> > I think, I should change something during minimization.
> > I've found this script from D. SIndhikara's tutorial:
> >
> > % minimization of solvent orientation
> >
> > &cntrl
> > ntmin = 1,
> > imin = 1,
> > maxcyc = 1000,
> > ncyc = 250,
> > ntb = 1,
> > igb = 0,
> > ntr=1,
> > restraint_wt=10.0,
> >
> > restraintmask='(:WAT.O) & ((!:WAT) <:5)',
> >
> > ibelly=1,
> > bellymask=':WAT',
> >
>
> This looks weird to me. Typically you use *either* restraints *or* belly,
> but not both (and I would suggest using restraints instead of belly).
>
> HTH,
> Jason
>
> --
> Jason M. Swails
> BioMaPS,
> Rutgers University
> Postdoctoral Researcher
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Oct 22 2015 - 07:00:07 PDT
