[Bug] How to set `pos0` for a harmonic trap?

[RodenLuo]

Describe the question A simulation can run fine without external forces. After adding only one harmonic trap to the first NT, using the center of mass coordinates from the configuration as pos0, the simulation explodes. How can I set the pos0 of a NT to be its initial position?

To Reproduce Run the simulation in pos0ForHarmonicTrap.zip.

Expected behavior The simulation runs successfully with a harmonic trap set on the first NT.

Desktop (please complete the following information):

• OS: Ubuntu 18.04.6 LTS
• oxDNA Version: RELEASE: v3.5.2 (GIT COMMIT: bf8729eb)
  • cmake/make/compiler versions: cmake version 3.21.1, gcc (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0

Additional context The structure comes from https://nanobase.org/structure/133. Thanks!

[ErikPoppleton]

Your force file looks more like a string than a harmonic trap. Are you trying to pull on your molecule or hold it in place? Also note that 1 is a very strong force stiffness, you might want to start with a much lower value and not have it increase over time (depending on what you're trying to model).

[RodenLuo]

Ah, very sorry. I attached the wrong force.txt file (probably some old force.txt in my Download folder). The one for the problematic run is here:

{
type = trap
particle = 0
pos0 = -40.6, -0.5, -5.6
stiff = 0.1
rate = 0
dir = 0.,0.,1.
}

#{
#type = trap
#particle = 7248
#pos0 = -43.31, 2.24, -0.47
#stiff = 0.1
#rate = 0
#dir = 0.,0.,-1.
#}

I am trying to reproduce the "harmonic traps" here, https://pubs.acs.org/doi/suppl/10.1021/acsnano.8b01844/suppl_file/nn8b01844_si_003.pdf (page 5, last paragraph). According to my calculations, using the units found here, page 13, I should set stiff = 0.1 and rate = 6e-6. Kindly correctly me if wrong: 1 unit of force constant = 57.09 pN/nm, the paper used 5.7 pN/nm, so the stiff = 0.1. 1 unit of length = 0.8518 nm and 1 unit of time = 3.03 ps gives 0.8518/3.03 = 0.28e12 nm/s. I'm using dt = 0.003. The paper's rate is 5.6e8 nm s−1. So I should set rate as 0.003 * 5.6e8 / 0.28e12 = 6e-6. I am a little bit inconfident about the handeling of dt = 0.003 as the documentation only says length simulation units/time steps.

In the force txt I give here, the rate = 0 because I see explosion when using rate = 6e-6 and then I am trying to debug it. I soon realize that even with a 0 rate, as long as I set the force, the simulation explode.

[RodenLuo]

The zip file in the initial post has been corrected.

[RodenLuo]

It is caused by the periodic bounding box. I extract the corrected position and make the simulation run. I have not figured out how to solve it gracefully yet. The information below might be helpful for further debugging.

Adding the following lines

                printf("ppos.x, ppos.y, ppos.z: %f, %f, %f\n", ppos.x, ppos.y, ppos.z);
                printf("tx, ty, tz: %f, %f, %f\n", tx, ty, tz);
                printf("ppos.x - tx, ppos.y - ty, ppos.z - tz: %f, %f, %f\n", ppos.x - tx, ppos.y - ty, ppos.z - tz);

in between Line 151 and 152 of the following code snippet

https://github.com/lorenzo-rovigatti/oxDNA/blob/a6214b99b64806d322e110c3d84fce4bbe9815f2/src/CUDA/Backends/CUDA_MD.cuh#L144-L153

gives the following output

ppos.x, ppos.y, ppos.z: -40.613007, 437.494629, 432.335693
tx, ty, tz: -40.599998, -0.500000, -5.600000
ppos.x - tx, ppos.y - ty, ppos.z - tz: -0.013008, 437.994629, 437.935699

Given pos0 = -40.6, -0.5, -5.6 and the first NT's position is -40.613006591796875 -0.50537109375 -5.664306640625, and the box sides are b = 438 438 438. It is clear that the simulator for some reason corrected the Y and Z of the first NT by one box side.

After setting pos0 = -40.6, 437.5, 432.3, the simulation runs fine.

[RodenLuo]

In the pos0ForHarmonicTrap.zip I attached, the box sides are b = 141 141 141. This is automatically shrinked by oxView when loading the .oxview file. I later changed it to the original value (438) in the .oxview file during debugging.

[lorenzo-rovigatti]

Thanks for the detailed info Roden. This is an issue that is halfway between a bug and a feature. The problem is that by default we internally always keep nucleotide coordinates within the box in order to improve numerical precision. However, while on CPU force computation is done on particles' absolute positions, on GPU this would decrease performance, so we always apply forces to in-box coordinates.

First of all, a workaround is to set fix_diffusion = false. Secondly, I'm thinking how to help users understand what's going on when they use forces that are computed with respect to positions that are outside of the box. I think the best option would be to not run any simulation, telling the user to either use in-box force positions, set fix_diffusion = false or CUDA_skip_force_position_check = false (which would be a new option). What do you guys think?

[RodenLuo]

on CPU force computation is done on particles' absolute positions

Does this mean on CPU, the default is fix_diffusion = false?

on GPU this would decrease performance

Does this mean, on GPU, if fix_diffusion = false, the performance will be decreased? If yes, is the decrease in terms of speed or in terms of precision?

I think I understand your proposed option and I believe it is a good solution.

What I still don't understand is, when the input configuration file has the first NT's position as -40.599998, -0.500000, -5.600000 (and the rest NTs' position are nearby), why the simulator believe the in-box-coordicates for the first NT is -40.613007, 437.494629, 432.335693. Why doesn't the simulator use the box that contain the original position data in the input configuration file?

[lorenzo-rovigatti]

Does this mean on CPU, the default is fix_diffusion = false?

No, by default fix_diffusion is always set to true. However, on CPU we have an additional data structure that keeps track of the times each particle crosses the box boundaries, which makes it possible for forces to act on absolute positions.

Does this mean, on GPU, if fix_diffusion = false, the performance will be decreased? If yes, is the decrease in terms of speed or in terms of precision?

The latter. On GPU we the most computational-intensive calculations are carried out with 32bit floating point arithmetic, whose numerical performance decreases as the absolute value of the numbers involved in the computations increases.

What I still don't understand is, when the input configuration file has the first NT's position as -40.599998, -0.500000, -5.600000 (and the rest NTs' position are nearby), why the simulator believe the in-box-coordicates for the first NT is -40.613007, 437.494629, 432.335693. Why doesn't the simulator use the box that contain the original position data in the input configuration file?

The box where particle coordinates are folded is always the "main" one. In other words, every time fix_diffusion is applied, each coordinate k of the centres of mass of each strand range between 0 and L_k.

[RodenLuo]

Thanks very much for the detailed explanation! I will leave this issue open. But please feel free to close it anytime from now till the potential code or documentation changes.

[RodenLuo]

I have realized that the "main" box can change during the simulation. So, it is impossible to "use in-box force positions" and set fix_diffusion = true, as one cannot set a fixed force position to guarantee it is always in the "main" box on the fly.

After re-reading your proposed solution, shouldn't the following line

CUDA_skip_force_position_check = false

be

CUDA_skip_force_position_check = true

or

CUDA_force_position_check = false?

I think a better solution would be to let the force positions be diffused the same way as the NT positions.

[lorenzo-rovigatti]

The main box can change if one simulates in the NPT ensemble, which is pretty unheard of for DNA origami simulations (as far as I know).

Yes, you are right, my post should read CUDA_skip_force_position_check = true, thanks for noticing!

[RodenLuo]

With the above attached zip file, change the force.txt to be the following so that the simulation would not explode in the very beginning.

{
type = trap
particle = 0
pos0 = -40.6, 437.5, 432.3
stiff = 0.1
rate = 6e-6
dir = 0.,0.,1.
}

{
type = trap
particle = 7248
pos0 = -43.3, 440.2, 437.5
stiff = 0.1
rate = 6e-6
dir = 0.,0.,-1.
}

The simulation explodes at t = 303.0000, which is simulation step 101000 as the dt = 0.003. The stiff at this step should be 0.606.

      294.0000  -1.508537   0.293211  -1.215326 
      297.0000  -1.510813   0.295916  -1.214897 
      300.0000  -1.510678   0.295441  -1.215237 
      303.0000  138188349.214137   2.919570  138188352.133707 
      306.0000  621847576.743570  63810.080298  621911386.823868 
      309.0000  22041041939.787567  485500552191059934838784.000000  485500552191081946546176.000000 
      312.0000  899053409797.718018  6323904809091587122030299840512.000000  6323904809091587122030299840512.000000 

The simulation settings have ensemble = nvt and do not include anything about fix_diffusion.

My guess is that after t = 300.0000, the main box changed.

[lorenzo-rovigatti]

ensemble is a key that is used only when simulating with MC, so it's not relevant here. If you haven't set fix_diffusion = false then you are experiencing the same exact issue as before.

[RodenLuo]

So, I was thinking about not setting fix_diffusion = false so as to not decrease the precision. The simulation runs fine for the first 99000 steps. And the first NT's pos is always near -44.239368, 439.182251, 436.191742, which is always near the harmonic trap I set. But then, at step 100000, the first NT's pos becomes 393.287781, 1.518220, 436.454224. See the below log for details. The box sides are b = 438 438 438. It apparently shifts to the nearby box. Is it possible to shift the traps pos at the same time? Is this a good idea or it does not make sense at all?

If you can kindly pinpoint the file and lines that detect and perform the fix_diffusion functionality, I might be able to implement the shift for the harmonic trap.

step: 99000
ppos.x, ppos.y, ppos.z: -44.239368, 439.182251, 436.191742
tx, ty, tz: -43.299999, 440.200012, 437.500000
ppos.x - tx, ppos.y - ty, ppos.z - tz: -0.939369, -1.017761, -1.308258
step: 100000
ppos.x, ppos.y, ppos.z: 393.287781, 1.518220, 436.454224
tx, ty, tz: -43.299999, 440.200012, 437.500000
ppos.x - tx, ppos.y - ty, ppos.z - tz: 436.587769, -438.681793, -1.045776

[lorenzo-rovigatti]

By default positions are wrapped into the main box size every fix_diffusion_every steps, which by default is 10^5. The fix_diffusion procedure is carried out in managers/SimManager.cpp, line 166.

[lorenzo-rovigatti]

Hey Roden, I commited a change that should take care of this issue. It is a bit more convoluted than yours, but it should cover more cases. It seems to work on your folder. Let me know how it goes, and feel free to reopen the issue if you find something wrong. I'll leave your PR open for now, in case my fix doesn't work.

[RodenLuo]

It runs through! Thanks for the fixing!