GROMACS version: 2019.4
Hi users,
I am trying to run a coarse grained simulation with Trp-cage (20 residues) in a modified CG water model box (the BMW model). The vdW interaction between the waters is not of LJ form so the user defined tables have to be provided. After solvation of Trp-cage in the water box, I ran an energy minimization. Up untill here the simulation looked fine. I checked the simulation box and there did not seem to have obvious overlapping or staggering. I then ran an NPT equilibration with this command:
“mpirun -np 27 gmx mdrun -npme 9 -noappend -s npt1.tpr -cpi state.cpt -tableb …/…/table_a/table_a*.xvg …/…/table_d*.xvg -tablep table_SOL_NON.xvg”
where I provided tables for angular potential, dihedral potential and the pair interaction. I’m pretty sure there is no issue with the angular potential and the dihedral potential table as I have used them many times and everything was fine. My simulation stopped at the very beginning (less than 5 s wall time), and the error message looks like this:
"Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.
-
Box (3x3):* -
Box[ 0]={Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.* -
Box (3x3):* -
Box[ 0]={ -nan, -nan, -nan}* -
Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.*
-
Box (3x3):* -
Box[ 0]={ -nan, -nan, -nan}* -
Box[ 1]={ -nan, -nan, -nan}* -
Box[ 2]={ -nan, -nan, Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.* -
Box (3x3):* -
Box[ 0]={ -nan, -nan, -nan}* -
Box[ 1]={ -nan, -nan, -nan}* -
Box[ 2]={ -nan, -nan, -nan, -nan}* -
Box[ 1]={ -nan, -nan, -nan}* -
Box[ 2]={ -nan, -nan, -nan}* -
Can not fix pbc.* -
Box[ 1]={ -nan, -nan, -nan}* -
Box[ 2]={ -nan, -nan, -nan}* -
Can not fix pbc.* -
-nan}* -
Can not fix pbc.* -
-nan, -nan}* -
Can not fix pbc."*
I have looked at the tpr file with gmx dump and saw that my initial box size was fine. I have tried changing tau_p to 5.0 but it still didn’t work. Currently I’m not sure where this issue comes from. Please let me know if you have any suggestion or comment. Thank you so much in advance!
The mdp file (It seems that I can not attach a file here):
*; VARIOUS PREPROCESSING OPTIONS = *
title = BMW-Martini
*; RUN CONTROL PARAMETERS = *
integrator = md
*; start time and timestep in ps = *
tinit = 0.0
dt = 0.020
nsteps = 5000
*; number of steps for center of mass motion removal = *
nstcomm = 1
*;comm-grps = DPPC WD *
*; OUTPUT CONTROL OPTIONS = *
*; Output frequency for coords (x), velocities (v) and forces (f) = *
nstxout = 0
nstvout = 0
nstfout = 0
*; Output frequency for energies to log file and energy file = *
nstlog = 1000
nstenergy = 1000
*; Output frequency and precision for xtc file = *
nstxtcout = 1000
xtc_precision = 100
*; This selects the subset of atoms for the xtc file. You can = *
*; select multiple groups. By default all atoms will be written. = *
*xtc-grps = *
; Selection of energy groups =
;energygrps = SOL NON
;energygrp_table = SOL SOL NON NON SOL NON
; please try not to change the following parts, they define the
; force-field interactions
; Force field parameters begin
cutoff-scheme = group
*; NEIGHBORSEARCHING PARAMETERS = *
*; nblist update frequency = *
nstlist = 10
*; ns algorithm (simple or grid) = *
ns_type = grid
*; Periodic boundary conditions: xyz or none = *
pbc = xyz
*; nblist cut-off = *
rlist = 1.4
; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = pme
*fourierspacing = 0.2 *
pme_order = 6
;rcoulomb_switch = 0.0
rcoulomb = 1.4
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 1.3
*; Method for doing Van der Waals = *
*vdw_type = User *
*; cut-off lengths = *
;rvdw_switch = 1.0
rvdw = 1.4
*; Apply long range dispersion corrections for Energy and Pressure = *
DispCorr = No
; Force field parameters end
*; OPTIONS FOR WEAK COUPLING ALGORITHMS = *
*; Temperature coupling = *
tcoupl = Berendsen
*; Groups to couple separately = *
tc-grps = System
*; Time constant (ps) and reference temperature (K) = *
*tau_t = 0.2 *
*ref_t = 300 *
*; Pressure coupling = *
*Pcoupl = Berendsen *
Pcoupltype = isotropic
*; Time constant (ps), compressibility (1/bar) and reference P (bar) = *
*tau_p = 5.0 *
*compressibility = 3e-5 *
*ref_p = 1.0 *
*; GENERATE VELOCITIES FOR STARTUP RUN = *
gen_vel = yes
gen_temp = 105
gen_seed = 473529
*; OPTIONS FOR BONDS = *
*constraints = none *
*; Type of constraint algorithm = *
constraint_algorithm = Lincs
*; Do not constrain the start configuration = *
unconstrained_start = no
*; Highest order in the expansion of the constraint coupling matrix = *
lincs_order = 4
*; Lincs will write a warning to the stderr if in one step a bond = *
*; rotates over more degrees than = *
lincs_warnangle = 30
Regards,
Peiyin