AW: Solvent box fluctuations

From: Götz, Alexander (alexander.goetz_at_mytum.de)
Date: Wed Mar 29 2017 - 03:37:59 CDT

Hi Josh,

thanks for your reply. I rotated the system to be aligned with the z-axis and activate the constantRatio option. The good news is, that it helped for the problem of single side thinning. The bad news is, that it didn't solved my problem. The box is still thinning (now on both sides). However, I observed this to be correlated with a rotating/tilting peptide and my currently chosen 2 kcal/mol seems to be significantly to small.

The peptide often reaches angles >45°. I hope this is not an issue with a wrong reference quaternion (I use (1, 0, 0, 0)). for the first, I have increased the the force constant to 5 kcal/mol for the rotation and reduced it to 1 kcal/mol for the translation. I am not sure if one of my problems is, that I only use a subset of the peptide for the rotation and translation constrain (only heavy atoms of the core region, so for a 31 residue peptide the 24 residues) but I will see how this works out. I also use restarting of the colvars now which I didn't before.... However, this requires some quick and dirty hand editing of the colvar.state file between equilibration and production. Is there any better/cleaner option for that? Otherwise I will solve it with a short python script in the future but maybe there is something in NAMD. I am not a fan of such "hacking" approaches. Just setting the other constrains force constants to 0 leads to significant performance decreases (from 110ns/day down to 75ns/day) because the other biasing forces are still computed. I also have activated wrapAll but I don't see any problem, because my constrained molecules is not translating between images.

I am also thinking about using the useconstantArea option for my x-y plane but I am not sure how this will affect the physics of my system.

Another thing I observed is that our old simulations which didn't use constrains and only 1fs integration didnt't used the vdWForceSwitching option, which activates the original CHARMM force switching (we use C36 as forcefield). We also activated this in our membrane simulations recently and had the same issues with box fluctuations, which could then finally be solve by the constantRatio switch.

I know a, these are a lot of things, but maybe you have some ideas on that?

Alex

________________________________
Von: Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov>
Gesendet: Montag, 27. März 2017 19:13:37
An: namd-l_at_ks.uiuc.edu; Götz, Alexander
Betreff: Re: namd-l: Solvent box fluctuations

Hi Alex,

Does your peptide need to be aligned along X? In membrane protein land, we use the "useConstantRatio" along with "useFlexibleCell" so that the growth/shrinkage of the X and Y axes are coupled together. If you aligned the long axis of your peptide along Z, using those options might eliminate the problems you have of the short dimensions changing independently and getting close contact between protein copies.

-Josh

On 03/27/2017 11:01 AM, Götz, Alexander wrote:
Dear NAMD users,

I am currently facing some troubles with my simulation and couldn’t find any more solution on the mailing list. So I hope you have a solution, idea of a good hint.

I want to simulate a 35 residue peptide in TFE/TIP3 (equal numbers) in a NPT ensemble. The box has an initial size of 104*60*60 Angstrom and the perfectly alpha-helical peptide is aligned along the X-Axis. We used this box in the past with peptides of up to 31 residues without any problems for several years. However, the new peptides will certainly unfold a little bite more than the old ones (approx. 6-8 residues) and therefore the 60 Angstrom of the y and z-axis will very likely become too short.

Originally the peptide is heavily tumbling in the box. To keep the box dimensions, I wanted to get rid of the rotational and translational movements of the peptide by using colvars. In the beginning I introduced a rotational and a translational restrains. As a result, I ended up with a heavily fluctuating/breathing box, what finnaly lead to the well-known partitioning error. Also with margin 2.0 this happens but what is much more problematic is that my box sometimes gets so small on one site, that there are only a few Angstroms between the peptide and the boundaries. This also stops me from just increasing the marign to values above 2.0.

Googling the mailing list, I found some tips for using slower rates of pressure coupling during equilibration. So I have set up the following equilibration protocol now (I only write down what I changed from the original CHARMM GUI protocols for NVT equilibration and NPT run in water solvents):

1. NVT, 295K, reassignFreq 500, reassingTemp 500, timestep 1.0, langevinDamping 5.0, run 500000

2. NPT 295K, FlexibleCell yes, wrapNearest off, margin 2.0, reassignFreq 500, reassingTemp 295, timestep 1.0, langevinDamping 5.0, langevinPistonPeriod 1000, lagevinPistonDecay 500, run 500000

3. NPT 295K, FlexibleCell yes, wrapNearest off, margin 2.0, reassignFreq 500, reassingTemp 295, timestep 1.0, langevinDamping 5.0, langevinPistonPeriod 750, lagevinPistonDecay 250, run 500000

4. NPT 295K, FlexibleCell yes, wrapNearest off, margin 2.0, reassignFreq 500, reassingTemp 295, timestep 1.0, langevinDamping 5.0, langevinPistonPeriod 200, lagevinPistonDecay 100, run 500000

5. NPT 295K, FlexibleCell yes, wrapNearest off, margin 2.0, reassignFreq 500, reassingTemp 295, timestep 2.0, langevinDamping 5.0, langevinPistonPeriod 100, lagevinPistonDecay 50, run 500000

The Production run works with the same settings like the last step of the equilibration, without reassigning the temperature every 500 steps. I commonly end up with a box fluctuation error after 5ns. I tried several parameters and also made a lot of changes in my equilibration protocol (It has been much worse in the beginning, so the protocol now is the best I could get until now) and also switched of the translational restrain because translation would not be a problem. For the equilibration run, I use the standard CHARMM GUI colvars on the backbone and side-chains together with my rotational ones. The rotational restrain is only on the backbone CA atoms of a core region of the peptide (approx. 24 residues in the center). For production I only stay with the rotational constrain (script shown in the end of this message). Changing the integration timestep from 2fs to 1fs didn’t made any difference. From what I can see in VMD (didn’t checked the absolute values), the box immediately starts breathing/oscillating after switching from NVT to NPT and the amplitude of the oscillation increases with time and piston decay and period. Both is not unexpected to me; however, I would expect the oscillation to slow down with ongoing time but I seems to increase.

I don’t have any more ideas what causes this problem. The worst case solution would be to increase the box size to an equally spaced one which would not require any constrain. However, this will significantly blow up (110k atoms vs 28k) and slow down my simulation which currently runs for approx. 110ns/day without fine tuning. I hope you have a tip or solution for me.

Best regards

Alex

#Peptide Translation and Rotation###

Colvarstrajfrequency 100
Colvarsrestartfrequency 100

### ROTATIONAL RESTRAINING ###
colvar {
  name rot

  orientation {
    atoms {
      atomsFile peptide.pdb
      atomsCol B
      atomsColValue 1.0
    }
    refPositionsFile peptide.pdb
    refPositionsCol B
    refPositionsColValue 1.0
  }
}

harmonic {
  colvars rot # Acting on colvar "rot"
  centers (1.0, 0.0, 0.0, 0.0) # Unit quaternion (no rotation)
  forceConstant 1.0
}

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:21:10 CST