Convert triclinic box from NAMD to LAMMPS - ERROR: Fix npt/nph has tilted box too far in one step - periodic cell is too far from equilibrium state

LAMMPS LAMMPS General Discussion
,
1

Dear LAMMPS users,

I have a system of a protein in a triclinic box that runs fine in NAMD. I have run equilibration for longer than 100ns with this system. However, when I convert the equilibrated structure to LAMMPS and try to run further equilibration I get the following error:

ERROR: Fix npt/nph has tilted box too far in one step - periodic cell is too far from equilibrium state (…/fix_nh.cpp:1214)

I have done the same before (equilibrated a structure in NAMD and converted its final trajectory to LAMMPS using charmm2lammps.pl) with a different structure but didn’t get this error message. This other structure was larger, but had an orthogonal box.

I’m stuck in this problem for several weeks and have already tried many different things. I hope someone here can give some advice on how to overcome it.

Some more information:

  • charmm2lammps exports a .data file with orthogonal box dimensions. However, I managed, following 8.2.3. Triclinic (non-orthogonal) simulation boxes — LAMMPS documentation, to compute LAMMPS box sizes and tilt factors that correspond to the triclinic simulation box of the converted pdb. I read the data file in vmd (topo readlammpsdata *.data) and the box looks ok.

  • I use box tilt large command before reading the data file (otherwise LAMMPS complains about the tilt factors).

  • I combine fix_lang and fix_nph

  • I tried to:
    reduce the timestep (from 2 to 0.1fs),
    increase minimization steps,
    reduce the dump frequency to 1 (so I can visualize the structure right after the first step, but LAMMPS gives error at the 1st step)
    increase the Pdamp of the nph

  • Change box from orthorombic to triclinic to reduce vacuum between periodic images
    at the end of this discussion Axel Kohlmeyer mentions that “some functionality does not work with triclinic boxes”. Am I hitting one of these functionality? Which?

Please, let me know if I can provide any other information to help people better understand the problem. I tried attaching a compressed file, but it says that new users can’t do that.

Thank you in advance for any comments/suggestion/help!

Amadeus

Beginning of data file:
LAMMPS data file. CGCMM Style. atom_style full. Created by charmm2lammps v1.9.2 on Wed Jan 18 06:07:58 PM EST 2023

  316058  atoms
  279667  bonds
  405516  angles
  632430  dihedrals
   33185  impropers
   15740  crossterms

      36  atom types
      65  bond types
     150  angle types
     359  dihedral types
      23  improper types

     -7.023765421836   32.762234578164 xlo xhi
 -0.809697875847   26.206920075237 ylo yhi
 330.857656154946   3047.891164180654 zlo zhi
  -7.3706424743813  -257.3278266819468  -28.2351575993900 xy xz yz

input file:
##########################################################################
############################# INITIALIZATION #############################
##########################################################################

units real
dimension 3
boundary p p p
atom_style full

##########################################################################
############################## FORCE FIELDS ##############################
##########################################################################

pair_style lj/charmmfsw/coul/long 12 14
bond_style harmonic
angle_style charmm
dihedral_style charmmfsw
improper_style harmonic

pair_modify mix arithmetic
special_bonds charmm

##########################################################################
############################### READ DATA ###############################
##########################################################################

box tilt large
fix cmap all cmap charmm36.cmap
fix_modify cmap energy yes
read_data sys_EQ_TER.data fix cmap crossterm CMAP
pair_coeff 29 31 0.154919 3.24019863787641 0.154919 3.24019863787641
#change_box all triclinic # xy final -7.2750966277473 xz final -237.6912825182827 yz final -26.8525803784517 remap
kspace_style pppm 1e-6 # requires coul/long pair_style

##########################################################################
########################## VARIABLE DEFINITIONS ##########################
##########################################################################

variable T equal 310 # Temperature of the system in K
variable P equal 1.0 # Pressure of the system in atm
variable dt equal 0.1 # Timestep in fs
variable EqNpt equal 50000 # No. of steps of NPT Equilibration
variable ThermoFreq equal 1000 # Freq. (in time steps) in which the thermodynamic info is printed on the screen
variable DumpFreq equal 1 #50000 # Freq. (in time steps) in which data is saved in dump files
variable ResFreq equal 50000 # Freq. (in time steps) in which restart files are generated
variable TempDamp equal 100*{dt} # LAMMPS manual fix npt variable PressDamp equal 1000*{dt} # LAMMPS manual fix npt

##########################################################################
########################## SIMULATION SETTINGS ###########################
##########################################################################

timestep ${dt}
neighbor 5.0 bin
neigh_modify every 1 delay 0 check yes

##########################################################################
############################# MINIMIZATION ##############################
##########################################################################

thermo 10
thermo_style custom step pe press fnorm fmax
thermo_modify flush yes

fix 1 all box/relax aniso 0.0
min_style cg
minimize 0.0 0.0 1000 10000 # 0.0 1.0e-10 1000 10000
unfix 1
write_restart restart.min1

###########################################################################
############################## EQUILIBRATION ##############################
##########################################################################

reset_timestep 0

thermo ${ThermoFreq}
thermo_style custom step ke pe lx ly lz vol temp press
thermo_modify flush yes

compute peratom all pe/atom

fix fix_H all shake 1e-8 100 0 t 1 2 3 4 5 6 7 8 9 10
velocity all create $T 1234567 mom yes rot yes

Langevin thermostat
fix fix_lang all langevin $T T {TempDamp} 48279
Nose-Hoover barostat
fix fix_nph all nph aniso $P P {PressDamp} nreset 1000

dump 1 all custom ${DumpFreq} dump_sys_EQ_eq1.lammpstrj id type x y z

restart ${ResFreq} sys_EQ_eq1.restart

run ${EqNpt}

2

Why not use an orthogonal box?

3

I agree with @srtee, if you can avoid using tilted box it is always better. If you must use a tilted box, then have a look at atomsk, this software is amazing at making a box less tilted.

Simon

4

The best course of action with this type of errors is to study the relevant part of the LAMMPS source code (indicated in the error message - (…/fix_nh.cpp:1214)). The code is as follows (line numbers are a little bit different as your LAMMPS is older):

The error is triggered if the absolute value of xy or xz becomes larger than the box length in x direction (multiplied by TILTMAX which is equal to 1.5) or when the absolute value of yz exceeds the box length in y direction (times TILTMAX) during a simulation with fix nph.

The solution is to rotate your entire box, so that the conditions above aren’t fulfilled.

2 Likes
5

I think the key here is “in one step”. The box will automatically be flipped if it doesn’t move too far in one step, see fix nvt command — LAMMPS documentation. I think you need to reduce your timestep or do some minimization with fix box/relax

The flip keyword allows the tilt factors for a triclinic box to exceed half the distance of the parallel box length, as discussed below. If the flip value is set to yes, the bound is enforced by flipping the box when it is exceeded. If the flip value is set to no, the tilt will continue to change without flipping.
1 Like
6

Thank you @mkanski and @stamoor for your suggestions/comments.

@mkanski
Indeed my xz values is larger than the box in x direction. I’m afraid that if I rotate the box I will reach one of the other conditions (xy or yz), but will see if I can do that.
The error indeed has something to do with the nph ensemble. I can run simulations in nvt.

@stamoor
I was already using fix box/relax when getting the error (I also reduced the time step from 2fs to 0.1fs and it didn’t help).

I will try some other things and post here when I have more to comment on that.
Thank you!

7

Hi @Amadeus,

as @srtee and @simongravelle before, I’ve seen your thread since a couple of days and do not understand the point in using a triclinic box for a protein in water. There might be, but without more detail on what you try to achieve here, it is impossible to tell. Also, not only is your xz tilt larger than your x dimension, as you stated, but so is your yz tilt with regard to your y dimension. It seems to me shady to attempt to get LAMMPS to run smoothly with such a geometry since it was designed to avoid it in the first place.

Since I see no reason not to use an orthorombic box, why not simply try change_box all xy final 0. xz final 0. yz final 0.? Without remapping your coordinates or changing the boundaries, there should be no overlapping problem. I am more familiar in making solids supercells that way but, as far as I can tell after a quick look at the documentation, the topology of your molecules should also be preserved. I would be happy to get a confirmation here to be honest, but quick tests confirmed me that image flags were updated.

You would basically have the same system but with shifted periodic conditions which should not change much with regard to the physics and dynamics (from a computational point of view, Axel Kohlmeyer’s answer in the thread you linked says it all).

1 Like
8

I know GROMACS lets you use triclinic cells to simulate (e.g.) truncated octahedra instead of cuboids to simulate fewer waters per protein.

I agree that in the context of protein in water, if a particular simulation result depended on triclinic vs orthogonal boxes, then the result is probably quite suspicious to me already, since nature is highly unlikely to pack proteins in such precise spatial repetition.

1 Like
9

Hi @Germain,

I understand your, @srtee’s and @simongravelle’s concern.
If nothing works I will end up moving to an orthorhombic/orthogonal box. I’ve done this before with a similar, yet different, system. However, I would like to keep this specific system triclinic (as I ran equilibration in NAMD).
Let me try to explain you why I want to keep the box triclinic:

This system is not a regular fully solvated protein. It tries to reproduce the collagen microfibrillar structure based on the 3HR2 PDB:
[1] # https://www.pnas.org/doi/10.1073/pnas.0502718103
RCSB PDB - 3HR2: Low resolution, molecular envelope structure of type I collagen in situ determined by fiber diffraction. Single type I collagen molecule, post rigid body refinement, 'relaxed'

Unlike most proteins, collagen is not found isolated and fully solvated in nature, and it does not completely fold to perform a specific function. The association of collagen molecules into fibrils/fibers provide collagen-based tissues some specific and unique properties.
When compressed in the specific triclinic UC given by [1], and periodically replicated in space through PBCs, the model reproduces the experimentally determined microfibrillar structure [1].
Since my model derives from the 3HR2 PDB structure I am trying to keep the triclinic unit cell so I can mantain the fibrillar structure.

Rotating the system does not change much with regard to the dynamics, but it may affect how the collagen molecules interact with their images. Maybe after a long equilibration there will be no visible difference between the triclinic and orthohombic model, but initially there is.

image
image2082×412 39.9 KB

image source: Tools for defect calculations — ASE documentation

Thank you all for your comments and suggestions!
I will insist a bit more with the triclinic box before moving to an orthohombic one though.

10

With proteins interacting with their images I am not sure if the pressure calculations will be correct – but that’s another problem.

LAMMPS is more restrictive (for good reasons) on the tilt factors allowed, so you need to “re-box” your simulation. You can reset the tilt factor by adding a suitable integer multiple of lx (= xhi-xlo) to xz so that xz falls within +/- 0.5 lx. Similarly ly and yz. The change_box command should be able to handle this.

For a more long-term solution, you should remap axes so that the longest axis is the x-axis (so z to x, x to y, and y to z). This maps xz and yz to xy and xz respectively, and since lx will now be very long these tilt factors will almost certainly be allowed. You will have to do this manually in the data file – the relevant equations are in “howto triclinic” in the documentation.

2 Likes
11

In these situations, I use Platoon to compute the primitive unit cell and to remap the atom positions. Be aware that sometimes the atom order is also shuffled, and a bit of scripting is required to fix that.

12

I once did progressive translational and rotational fitting of a (plain text) trajectory in awk.

I will never do that in awk again.

EDIT: I joke about never wanting to awk again, but if you’re doing something rote to every row of a tabular data file and you’ve got the “programming mindset”, it’s ridiculously powerful and fast and found almost everywhere, including on HPC clusters.

1 Like
13

Hi everyone,

issue fixed!
As pointed out by @mkanski, the error was triggered because the absolute value of the xz tilt factor was greater than 1.5 times the box length lx.
I remapped the atoms to a new coordinate system such that the former long axis, the c vector of the fractional/crystal coordinate system, is now aligned with the x axis. I wrote an Matlab script to compute the rotation matrix and a Tcl script to rotate/remap the atom coordinates in vmd.

The new simulation box has a very long lx dimension, small tilt factors, and runs without triggering the error I mentioned in the post.

Thank you all for your suggestions and comments!

4 Likes