Re: [lammps-users] Conversion check from Gromacs to Lammps

LAMMPS LAMMPS Mailing List Mirror
1

As Axel said, it’s probably not a good idea to use a large complicated system to test your conversion. It’s better to use a tiny molecule for comparison. Create a simulation that contains only a single small molecule. (In the example below, I use ethylene.) Make sure this small molecule contains both dihedral and improper interactions. Make sure not to use a conformation of the molecule which is not in equilibrium. (In an equilibrium conformation, the forces on all the atoms will be 0. Probably both LAMMPS and GROMACS will agree that the forces are 0. This does not tell you anything.)

Question for Xuejiao:

Did you convert the entire GROMOS force field into moltemplate (LT) format?
Or did you convert only the ITP files containing the atomic coordinates and bond topology of a molecule into moltemplate format? (If you have code which can convert ITP files into moltemplate format, please share it with me, and I will find a way to give you credit for the contribution. I recently submitted a paper on moltemplate with all of the people who contributed force fields as coauthors.)

Andrew
P.S. If it helps, here is the protocol we used to test the DREIDING force field using ethylene. This is the way we did it, but I worry the long explanation below may be more confusing than helpful. Feel free to skip the rest of this email.

----- Details (feel free to skip) -----
When Matt Bone converted the DREIDING force field into moltemplate format, we ran into the same problem. We created some moltemplate files for ethylene in different conformations. If it helps, I have attached these files. Ethylene is probably the simplest and smallest molecule that you can use for this purpose. (Don’t use alkanes because they lack improper interactions.)
To test whether the dihedral interactions were converted successfully, calculate the forces on a conformation of ethylene which has been twisted around its central axis using both LAMMPS and GROMACS. (Although potential energy is not a meaningful way to compare the behavior of the molecule in LAMMPS and GROMACS, differences in energy are meaningful and can be compared. So if you don’t like measuring forces, you can instead measure the difference in potential energy of the twisted conformation of ethylene with the relaxed conformation of ethylene, using both LAMMPS and GROMACS. These energy differences should agree. The relevant files are “ethylene_rotate.lt” and “ethylene.lt”)
Similarly, to test the improper interactions, I would use a conformation where the 3 atoms surrounding each carbon atom are not coplanar. (Such a conformation is stored in the “ethylene_improper_sym.lt” file which is attached. I would try this first.)

Again, the moltemplate files I attached contain different conformations of the same ethylene molecules, but these were intended to be used with the DREIDING force field. You are not using the same force field, so you will have to create your own version of the ethylene LT file and copy the coordinates from my attached files into your LT file. You will also have to create equivalent versions of ethylene in GROMACS format so that you can compare with gromacs.

Protocol:

  1. Edit the “ethylene.lt”, “ethylene_improper_sym.lt”, and “ethylene_rotate.lt” files to make sure that the partial charge of each atom agrees with the charges that GROMACS is using. (Alternatively, set all atom charges to 0.)

  2. Edit the “system.lt” file and select which conformation of the molecule you want to simulate by uncommenting the line containing the version of the molecule you want. (Either “ethylene_improper_sym.lt”, or “ethylene_rotate.lt”, or “ethylene.lt”. I would start with “ethylene_improper_sym.lt”.)

  3. Run moltemplate:
    moltemplate system.lt

  4. Run LAMMPS
    lmp_mpi -i run_print_pe.in

(I created the “run_print_pe.in” script to make sure you are measuring the potential energy at t=0, before any relaxation has occurred. Note: The name of your LAMMPS executable binary (“lmp_mpi”) may vary.)

  1. Make note of the potential energy of the system when the molecule is not flat:

  2. Edit the “system.lt” file again, and uncomment the line containing “ethylene.lt”, and comment out the other lines. (ie, comment out the line containing “ethylene_improper_sym.lt”.)
    Then repeat steps 3-5)

(In the “ethylene.lt” file, the conformation of the molecule is flat and very close to its equilibrium conformation.)

  1. Subtract the potential energies (flat vs. not-flat).

  2. Do the same thing in GROMACS

README_visualize.txt (2.85 KB)

ethylene_improper_sym.lt (757 Bytes)

ethylene_rotate.lt (764 Bytes)

ethylene.lt (764 Bytes)

system.lt (211 Bytes)

run_print_pe.in (958 Bytes)

2

I worked with Gromacs initially but now I want to transfer to LAMMPS to use some of LAMMPS’ built-in functionalities. I therefore need to convert force field parameters and topology from Gromacs to LAMMPS. I am using LAMMPS-5Jun19. I did the conversion from Gromacs to LAMMPS with the help of moltemplate. Now I need to check whether the conversion went well or not. I compared the total potential energy and the individual parts of potential energy for the same system under NPT ensemble between two softwares. The total potential energy and individual parts were not completely identical, so I am writing this email to ask whether it is necessary to make sure these values to be identical or is there an acceptable difference (for example energy difference <<1%)?

that depends on many factors. differences can come from multiple sources:

  • floating point math is not associative: the order of floating-point math operations can impact the results. that is more apparent when the result of the some would be a small value due to summing terms with different signs (like when summing forces).
  • differences in the implementation of force kernels: there may be differences in whether and how a smoothing function is applied to forces around the cutoff
  • differences in the implementation of the neighbor list builds: there may be differences in how strict a cutoff is applied when building neighbor lists
  • differences in applying cutoffs: some MD codes apply cutoffs based on charge groups (= all atoms inside such a group contribute to forces and energy, even if some of them would be outside the cutoff) others use a strict cutoff.
  • approximations to expensive calculations: MD codes may use different approximations to speed up time critical parts of the calculation. this could be approximations to square roots, the complementary error function, or using mixing single precision with double precision floating point operation for better vectorization.

please note that gromacs is often built using single precision floating point math while LAMMPS uses double precision, gromacs favors speed and efficiency while LAMMPS favors flexibility and being able to handle very large systems. that has consequences for both, simulation speed and precision of the results.

I checked the mailing list for similar questions and I noticed that Axel mentioned more than once that it is meaningless to compare the absolute energies. I am wondering if that is because the potential energy is dependent on the positions of the atoms in the system? If yes, the potential energy at 0 timestep of the same system should be comparable, is that correct? If not, what other possible reasons make the potential energy incomparable?

total energies are not comparable because they can be computed differently. constant terms may be included or not without affecting the forces. for example, the potential energy may be shifted so that the energy at the cutoff would be zero (or not). it can also differ in the implementation of long-range electrostatics. again, constant terms may be included or not without affecting forces.

Axel also suggested to compare forces. The force per atom and the total force of the system at a specific timestep can be obtained for both softwares. For my system, the force per atom varied a lot from one software to another and the total force of the system was close to 0 for both softwares (but the values were not identical). It is however difficult for me to figure out what caused the big difference for the force per atom since I found out that the individual force parts (such as VdW force, electrostatics force, bond force and so on) are not easy to be obtained as the individual potential energies. I am wondering if I have to compare the forces to confirm my conversion is good, what are the corresponding requirements? Is the total force of the system close to each other is enough? Or the force per atom difference indicates that the conversion is not good? If yes, how can I track back to the causes? Or maybe someone has a better idea to help me to check whether the conversion is good or not?

this is all easy to compare and test by not testing your “production” system, but by constructing test examples where you can deliberately turn of contributions to the potential functions, e.g. by (selectively) setting the LJ epsilon or individual charges to zero. and again, keep in mind that when comparing results of a single precision gromacs binary to a double precision LAMMPS binary, there are differences to be expected and they are expected to grow exponentially. so the best is to compare computed properties at step 0.
axel.

As Axel said, it’s probably not a good idea to use a large complicated system to test your conversion. It’s better to use a tiny molecule for comparison. Create a simulation that contains only a single small molecule. (In the example below, I use ethylene.) Make sure this small molecule contains both dihedral and improper interactions. Make sure not to use a conformation of the molecule which is not in equilibrium.

Typo.
I meant to say “Make sure to use a conformation of the molecule which is not in equilibrium.”