hoomd.external.field.Periodic in 2-dimensions with Cell type neighbor list

85 views
Skip to first unread message

kmo...@ucsb.edu

unread,
Jan 28, 2023, 4:30:04 PM1/28/23
to hoomd-users
Dear All,

I am attempting to run 2D simulations with Brownian particles under an external periodic field. I believe that there is an error in the implementation of the external periodic potential in 2D, at least with the cell type neighbor list. As soon as I implement a potential, no matter how weak, there is an error computing the neighborlist and my particles have position NaN.

I've attached a small simple script to reproduce the error. I've noticed that there is no error when I switch  to the Tree type neighborlist, but unfortunately the system size I am interested in runs faster with the Cell type.This problem also does not appear in 3D simulations.

Is this some subtle problem with my implementation, or an error with the code. I am running Hoomd v3.8.0

Thank you for any advice you can give,
Kevin

examplePeriodicFieldError.py

Joshua Anderson

unread,
Jan 30, 2023, 8:29:47 AM1/30/23
to hoomd...@googlegroups.com
Kevin,

Yes, `V_box` will be 0 in 2D simulations:
which will result in a divide by 0:

You should be able to confirm the presence of `NaN` forces in `periodicForce.forces` after a call to run(0). No other forces, neighbor lists, or integrators are needed in a minimal script that reproduces the behavior you are seeing. Look closer at the 2nd code snippet (L165-169) and you will see that the computation involves the cross product of the 3D box vectors. To fix the behavior you describe, you can either (or both):

a) File a bug https://github.com/glotzerlab/hoomd-blue/issues/new?assignees=&labels=bug&template=bug_report.yml regarding the NaN forces. The fix would be to detect 2D boxes and raise an exception instead of computing NaNs.

b) File a feature request asking for 2D support in `md.external.Periodic` https://github.com/glotzerlab/hoomd-blue/issues/new?assignees=&labels=enhancement&template=feature_request.yml. As the first and user asking about this, you would be the one making the changes to the code and submitting the pull request. See the contribution guidelines here https://hoomd-blue.readthedocs.io/en/v3.8.1/contributing.html and we are happy to answer specific questions on the GitHub issue comments thread(s).
------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/e7f372ba-86a7-4585-8910-52b3682b1301n%40googlegroups.com.
<examplePeriodicFieldError.py>

kmo...@ucsb.edu

unread,
Jan 30, 2023, 2:53:24 PM1/30/23
to hoomd-users
Dear Joshua Anderson,

Thank you for explaining the problem. I apologize if this is a naive question, I'm not particularly skilled at reading C++ code, but shouldn't the divide by zero be taken care of in the BoxDim.h file?

I see there is a flag "twod" that causes the method "getVolume" to return the area if true. I'm not exactly sure how the parameter twod gets tracked in the evaluator, but comparing with the older code v2.9 when it used to work, the actual evaluation of the force energy and virial are practically identical.


I'll work on submitting a formal bug or feature request shortly, but I just have to parse through some C++ code to understand exactly what I'd be changing first.


Thank you,
Kevin

Joshua Anderson

unread,
Feb 3, 2023, 8:41:48 AM2/3/23
to hoomd...@googlegroups.com
Kevin,

BoxDim does not know whether the user is simulating a 2D or 3D system, so the caller to getVolume must provide the information. Yes, it is a poor design, but it would be take a large amount of effort to fix across the whole code base.

I don't see how external.periodic code could have worked in v2. In v2, users were encouraged - but not required - to set Lz=1 in 2D boxes. BoxDim computes the volume as Lx * Ly * Lz, so as long as Lz = 1, the value given is the correct numerical value for the box area. However, the remainder of the external.periodic code computes cross products and operates on the box as if it has 3D periodicity. At first glance, I don't see how this could be correct in 2D. One would need to write a completely separate code path.

When using v3, Lz will be 0 in most cases for 2D boxes. However, this is not a given. There are some corner cases (to support reading gsd files written by HOOMD v2) where a non-zero Lz is accepted, but ignored for 2D simulations in v3.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
Reply all
Reply to author
Forward
0 new messages