Strange electrostatic potential THF

[al4869]

GROMACS version: 2020.5 & 2021.1
GROMACS modification: No

Hello,
I am working on a box (20x20x20 Angstrom) containing several (in line with the experimental density) tetrahydrofurane (THF) molecules randomly arranged in the space (neutral molecules), and, after performing NVT MD on such a system (more details below), I would like to compute how the electrostatic potential changes along z-axis.
For this, I used

 gmx potential -s NVT_md.tpr  -f NVT_md.xtc  -o potential.xvg -sl  10000

Now, since the molecules are randomly arranged, I would expect a rather fluctuating potential like the one in panel A below (that I obtained for a NVT MD on a box composed by pentacene molecules)

Quite on the contrary, I actually obtain a monotonically increasing potential curve for the THF system, like the one in panel B above

I have tried to understand possible causes for this strange behaviour in the case of THF but:

• Calculating the potential on a snapshot or averaging over up to 50000 steps of the trajectory leads to qualitatively similar results
• Calculating the potential on wrapped or unwrapped box (since PBC are not taken into account in “gmx potential” as per the manual) leads to qualitatively similar results
• Using different force field (OPLS or even Quantum-Mechanically derived force fields) leads to qualitatively similar results
• Performing NPT or NVT MD leads to qualitatively similar results
• Using larger or smaller boxes leads to qualitatively similar results
• Using different versions (I tried 2020.5 and 2021.2) of Gromacs leads to qualitatively similar results
• Finally, I have also tried to perform MD adding vacuum either (i) on the right of the z-axis of the box or (ii) both on the left and the right of the z-axis of the box (and imposing constraints to the molecules along z direction so that they do not move inside the “vacuum” region), but the results are even stranger:

The image in panel C refers to the case where the vacuum is only on the right of z-axis: you can see that the potential fluctuates and increases in the region with the molecules and then increases as a straight line in the vacuum region.

The image in panel D below refers instead to the case with vacuum both on the left and on the right of z-axis (the molecules are in the center). Once again you can see that the potential start increasing in the region (2.5 to 4.5 nm) where the molecules are located, but strangely this time there is no increase in the vacuum region.

I have also visually inspected the snapshot and it seems that there are no strange arrangements of the molecules (moreover, I got similar results also analyzing directly the structure obtained arranging THF molecules through PACKMOL, i.e. completely randomly arranged)

I have browsed the forum and also the web, but I have not found anything similar to this result. I also asked to other colleagues (working in bioinformatics, so proteins in water) and they have never obtained something like this. It seems to be related to THF, but I honestly do not understand why: I am quite new to Gromacs, and probably I am missing something, but I cannot think to any other analysis to spot the problem, so I am asking for help to this forum where many people has much more experience than me.

I attach the files I used for the NVT simulations using OPLS (a zip folder, but addingextension .txt so I can upload it) if needed.

thf.zip.txt (20.7 KB)

Many thanks for any help you can give me,
Alessandro

[mick]

Out of curiosity, is there any particular reason why you are restricted to looking at the Z-axis? What do the X and Y look like?

Have you tried opening the snapshot in a different program (say VMD) and computing the potential there to compare?

The VMD PMEPot Plugin should work for this, but you’ll need to output a pdbqr file. You can do this using gmx editconf and writing a loop to replace your “q” column with the charge values listed in your topology file.

[hess]

You have not written how you treat the electrostatic interactions, that might have quite some effect.

Also you box is very small. I would suggest to use at least 4x4x4 nm and equilibrate sufficiently long.

[al4869]

Thanks for the reply. Along X and Y the values are different, but the trend is the same, as you can see in the following

I will try to use VMD plugin to compute the potential and let you know if this solves the issue

[al4869]

Thanks for the reply.
I attached to the first message all my .mdp files, to share my settings. However, here is the part concerning the electrostatics:

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = cut-off
rcoulomb-switch          = 0
rcoulomb                 = 0.7
; Dielectric constant (DC) for cut-off or DC of reaction field
epsilon_r                = 1
; Method for doing Van der Waals
vdw-type                 = Cut-off
; cut-off lengths
rvdw-switch              = 0
rvdw                     = 0.7
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = no
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
 fourier_nx               = 0
 fourier_ny               = 0
 fourier_nz               = 0
; EWALD/PME/PPPM parameters
pme_order                = 4
ewald_rtol               = 1e-05
epsilon_surface          = 0
optimize_fft             = no

Concerning the fact that the box is small, in my first post I wrote explicitly that “Using larger or smaller boxes leads to qualitatively similar results”. I have tried also 5x5x5 nm and the trend is similar

[hess]

You are using a plain cut-off for Coulomb with a cut-off distance of 0.7 nm. This will give extreme artifacts. Use PME instead.

[al4869]

Thanks for the reply.
If you mean replacing

@al4869

coulombtype = cut-off

with

coulombtype              = PME

I forgot to mention it in the original post, but I tried this either and the trend is once again similar:

[hess]

How long did you simulate? I wouldn’t be surprised if the instantaneous field looks as it does there. In the limit of infinite sampling it will be zero.

[al4869]

After a minimization, I equilibrated for 1 ns and then run a production MD of 1 ns, printing 50000 snapshot in the trajectory file. Do I need more than 50000 snapshots, in your opinion?

[hess]

No, that’s more than sufficient snapshots. But you might need much more than 1 ns.

[al4869]

50 ns, the problem persists

Plus, I have asked to a colleague and she found a similar trend for a 500-ns simulation in water
At this point I highly suspect a bug in the code

[al4869]

Here is the potential for the 500-ns simulation of a protein in water I was talking about in the previous post

[hess]

At 50 ns the potential is lower by nearly an order of magnitude.

A protein can have a large dipole and has an extremely long (on MD time scales) rotational diffusion time.

I don’t see any issue here.

[al4869]

Maybe I am missing something.

Below, I report the result of the potential from a DFT-MD over a THF snapshot with vacuum along z-axis. The potential oscillates in the part where molecules are present and then is almost flat in the “vacuum” part, as I would expect. If compared with panel C of my first post (same system, but with potential evaluated through gmx potential) I do see a strong difference in the trends.

Since I get similar behaviour (monotonous change in the potential over the box) also for almost apolar solvents (I simulated a box of iso-octane molecules), maybe I am missing something, so I will continue investigating.

[hess]

I assume there was a requirement of the potential being equal on both sides of the box in the DFT result. In gmx potential there is no such requirement. Adding that requirement will likely make the potential flat in the vacuum. Another aspect that has a large effect is the boundary condition at infinity for PME. Default is conducting boundary conditions, which results in a far larger dipole, and thus field, over the box. Using insulating boundary conditions, epsilon-surface=0, suppresses the net dipole by a factor of the dielectric constant.