Lennard Jones Shifted Radius

User discussions
1

GROMACS version: 2018
GROMACS modification: Yes/No

I am trying to shifted the Lennard Jones interaction - this means: instead of U(r), I want an interaction of U(r-ro) where ro = 0.2 nm (for example).

I have a code in order to generate the tables according to the GROMACS documentation. The table contains x,f(x),f´(x),g(x),g´(x),h(x),h´(x).

For ro = 0 nm and ro = 0.1 nm, the simulation is working. However, when I tried to run the case of ro = 0.2 nm, the minimization is blowing up.

This is the error:

Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 100 (which may not be possible for your system).
It stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.

Double precision normally gives you higher accuracy, but this is often not
needed for preparing to run molecular dynamics.
You might need to increase your constraint accuracy, or turn
off constraints altogether (set constraints = none in mdp file)

writing lowest energy coordinates.

Steepest Descents converged to machine precision in 22158 steps,
but did not reach the requested Fmax < 100.
Potential Energy = 1.1675758e+11
Maximum force = 4.6688879e+10 on atom 16626
Norm of force = 5.9937074e+08

The Fortran code that I am using is:

  *program gentable*

  •  implicit none*

  •  real, parameter :: delr=0.005, rcut=1.2, small_r_threshold=3.8e-2*

  •  real :: r, term1, term2, term3, term4, term5, term6, term7*

  •  integer :: nbins, j*

  •  nbins = int( (rcut+1)/delr ) + 1*

  •  *

  •  open(10, file='table_2.xvg', status='replace')*

  •  do j = 0, nbins*

  •      r = delr * j*

  •      if (r < small_r_threshold+0.2) then*

  •          ! If r is very small, set the terms to zero to avoid division by zero*

  •          term1 = 0.0*

  •          term2 = 0.0*

  •          term3 = 0.0*

  •          term4 = 0.0*

  •          term5 = 0.0*

  •          term6 = 0.0*

  •      else*

  •          ! Normal calculation for larger r*

  •          term1 = 1.0 / (r)*

  •          term2 = 1.0 / ((r)**2)*

  •          term3 = -1.0 / ((r-0.2)**6)*

  •          term4 = -6.0 / ((r-0.2)**7)*

  •          term5 = 1.0 / ((r-0.2)**12)*

  •          term6 = 12.0 / ((r-0.2)**13)*

  •      end if*

  •  write(10,*) r, term1, term2, term3, term4, term5, term6*

  •  end do*

  •  end program*

Any suggestions? I don´t understand why is working for 0 nm and 0.1 nm but nor for 0.2 nm

Thanks!

2

My guess is that there is at least one distance smaller than 0.2 nm which makes the system blow up. I would suggest to use a constant force, and thus a quadratic potential, for x < 0.2 + epsilon, where epsilon is a small number.

3

hey Hess! thank you for your reply. I can’t use a quadratic potential because I want to use the LJ potential interaction with a shift in the radius.

I don’t know if I have to consider something else (for this case) when I generated the table with the forces and the LJ potential because it worked in the previous two cases.

Also, if I increment the size of the box by 100 times, it is not blowing up - weird

4

What I wrote is that you should use a quadratic potential only at r < 0.2+eps, e.g r < 0.201. That you can do.

5

I tried it but it is blowing up

6

But how do start the simulation? If you shift the LJ potential by 0.2 and do not scale the starting coordinates you create an enormous pressure and the system will likely explode.

7

Sorry Hess - I been out of the loop! Yes, you are right. One question: How can I scale the starting coordinates?

8

Editconf can scale the coordinates. But if you have bonded interactions, bonded distances will also be scaled.

9

Hi hess. I am sure there is no overlapping atoms and the LJ shifted potential (custom table) is applied to the interaction between CH4 and SOL molecules. However, I think the minimization is blowing up because for the SOL SOL interactions is using this custom table – I want that to use the TIP4P force field for the SOL SOL interaction. Do you know how I implement this correctly? Or any suggestions? Thanks!

I check the energies after the minimization:

Screenshot 2025-03-14 at 4.56.04 AM
Screenshot 2025-03-14 at 4.56.04 AM1220×540 174 KB

PS: I am attaching my topol.top file and my system is a cubic box with 1 CH4 molecule in the center of the box and 6000+ SOL molecules.

topol.top (2.2 KB)

10

You will need to generate a standard LJ+Coulomb table for SOL-SOL.

11

perfect, I was thinking about it! thank you Hess!