Project

General

Profile

Bug #177

position and velocity set to nan

Added by Yang Ye about 12 years ago. Updated almost 12 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

it is an customized FF with section

[ dihedraltypes ]
; i l func q0 cq
;added by Dongsheng
OM OM 0 0.000 1 2
OA OA 0 0.000 1 2
OM OA 0 0.000 1 2

in ffpolymerbon.itp

If there is entries like following in topol.top
[ dihedrals ]
1 2 3 4 1

It will result in nan in position and velocity and cause error
-------------------------------------------------------
Program mdrun_s, VERSION 3.3.1
Source code file: nsgrid.c, line: 220

Fatal error:
Number of grid cells is zero. Probably the system and box collapsed.
-------------------------------------------------------

I have verified in bondfree.c that coordinates are the same after calling pdihs.

If there is no dihedral defined in topol.top. The simulation is alright.

Is there any other place which changes the coordinate?

temp2.tar.gz (24.1 KB) temp2.tar.gz test case Yang Ye, 11/16/2007 11:23 PM
vec.h (19.7 KB) vec.h gmx/include/vec.h David van der Spoel, 11/19/2007 07:19 PM
bondfree.c (47.1 KB) bondfree.c gmx/src/gmxlib/bondfree.c David van der Spoel, 11/19/2007 07:19 PM

History

#1 Updated by Yang Ye about 12 years ago

Created an attachment (id=257)
test case

test file including force field and run files.

type "grompp" and "mdrun" to run the case.

#2 Updated by David van der Spoel almost 12 years ago

This is because your dihedrals are exactly 0, that is all 4 atoms are co-linear. With more "physical" input this wouldn't happen.

#3 Updated by David van der Spoel almost 12 years ago

OK, I have made a fix where the dihedral gives zero force if it is completely linear. That means the simulation will not crash but the forces may be discontinuous close to a linear molecule. In addition the linear angles will also give you a discontinuous force close to 180, basically because we disallow 180 degrees angles altogether. This is not a pretty solution but it's the best I can do.

The dihedral angle is defined as the angle between two planes, however we cannot make these planes since the vectors making up the planes are on a line, hence the planes are undefined.

Erik, shall I commit this fix?

#4 Updated by Erik Lindahl almost 12 years ago

Hi,

That sounds like the only reasonably alternative, so please do!

#5 Updated by David van der Spoel almost 12 years ago

Created an attachment (id=265)
gmx/include/vec.h

New vec.h

#6 Updated by David van der Spoel almost 12 years ago

Created an attachment (id=266)
gmx/src/gmxlib/bondfree.c

#7 Updated by David van der Spoel almost 12 years ago

With the two file attached the problem is solved, but please keep in mind my previous comments.

#8 Updated by Yang Ye almost 12 years ago

Thanks. I am aware for the effects of the fix.

Also available in: Atom PDF