normal NVT MD works but not NPT MD (domain decomposition error)

  • payel
  • payel's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 4 months ago #1299 by payel
Hello, I am trying to run MD of a 42 residue peptide solvated in water with Martini forcefield. I used the latest martinize script with the options -p backbone -elastic -ef 500 -el 0.5 -eu 0.9 -ea 0 -ep 0 and generated the needed files. I run a short minimization followed by a backbone restrained NPT MD with no problem. Now for the production run, if I use NPT conditions, then I get the following error message:

The charge group starting at atom 2178 moved than the distance allowed by the domain decomposition (1.200000) in direction Z
distance out of cell nan
Old coordinates: 0.259 0.237 0.529
New coordinates: nan nan nan
Old cell boundaries in direction Z: 0.000 2.839
New cell boundaries in direction Z: 0.000 2.866

Program mdrun_mpi_bg, VERSION 4.5.4
Source code file: domdec.c, line: 4124

Fatal error:
A charge group moved too far between two domain decomposition steps
This usually means that your system is not well equilibrated.

This happens after a few hundred ps of simulation. But I have no ptoblem running a NVT simulation for 8 ns. If I try to use the end structure of NVT run as the starting structure of NPT, same error occurs. I tried smaller timestep (0.02 ps). Below is the mdp file.

VARIOUS PREPROCESSING OPTIONS =
title = Martini
cpp = /usr/bin/cpp

; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0.0
dt = 0.020
nsteps = 500000
; number of steps for center of mass motion removal =
nstcomm = 1
comm-grps = System

; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 5000
nstvout = 5000
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 1000
nstenergy = 100
; 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 = Protein W ION

; 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.2

; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Shift
rcoulomb_switch = 0.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 15
; Method for doing Van der Waals =
vdw_type = Shift
; cut-off lengths =
rvdw_switch = 0.9
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = Berendsen
; Groups to couple separately =
tc-grps = Protein W ION
; Time constant (ps) and reference temperature (K) =
tau_t = 1.0 1.0 1.0
ref_t = 300 300 300.0
; Pressure coupling =
Pcoupl = berendsen
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p = 1.0
compressibility = 5e-5
ref_p = 1.0

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel = no
gen_temp = 300
gen_seed = 473529

; OPTIONS FOR BONDS =
constraints = none
; Type of constraint algorithm =
constraint_algorithm = Lincs
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


I am using gromacs 4.5 in Bluegen/P (64 processor). The system contains 10735 atoms in 10.8x10.8x10.8 nm3 box. Can someone is the forum can please help me?

Thanks a lot in advance.

Please Log in or Create an account to join the conversation.

More
11 years 4 months ago #1300 by djurre
There is nothing obviously wrong with your mdp-file.

After how many steps did it happen? Did you retry the simulation with different velocities? (gen_vel = yes). Does it happen again?

Please Log in or Create an account to join the conversation.

  • payel
  • payel's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 4 months ago #1301 by payel
Thanks for your reply. Well, I guess I found a way out. The error occured around ~40000 steps when I used domain decomposition. But I just managed to run it upto ~175000 steps using particle decmposition. And then I get the typical Lincs warning error for a ILE residue.

Step 174877, time 5246.31 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms nan, max 136.105209 (between atoms 74 and 75)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
45 46 33.0 0.2650 0.2650 0.2650
72 73 90.0 0.3100 1.6004 0.3100
74 75 90.0 10.1973 42.5026 0.3100
81 82 90.0 6.5496 22.9088 0.2650
Wrote pdb files with previous and current coordinates


Can you suggest something? Should I change to long stiff bonds instead?

Thanks again.

Payel

Please Log in or Create an account to join the conversation.

More
11 years 4 months ago #1304 by djurre
I would try to use domain decomposition if possible, since it is much faster. If it doesn't work and keeps giving crashes that probably means something is wrong with the way you setup your system.

You are using Martini with a general purpose elastic network. Are you sure you don't want to use Elnedyn (one of the forcefields available in martinize.py)? That works well for most people. Also you said you have small peptide. Maybe you can do without elastic network all together?

Please Log in or Create an account to join the conversation.

  • payel
  • payel's Avatar Topic Author
  • Offline
  • Fresh Boarder
More
11 years 4 months ago #1306 by payel
Djurre,

Thanks for your reply. Changing constraints for BB/SC-SC interaction to stiff bonds fixed the problem. I am using elastic network to include secondary structure constraints. May be I will try elnedyn later.

Payel

Please Log in or Create an account to join the conversation.

More
11 years 4 months ago #1307 by djurre
Secondary structure is already conserved in Martini by means of angles and dihedral angles. You specify the secondary structure to martinize.py, either using "-dssp DSSP.exe" or giving a string of letters. For that you do not need an elastic network.

The elastic network you need for ternary structures.

Please Log in or Create an account to join the conversation.

Time to create page: 0.097 seconds