Project

General

Profile

Bug #590

Problem with PME and Test particle insertion

Added by Kevin Daly about 9 years ago. Updated about 9 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

Recently I ran some test particle insertion simulations of TIP4P/2005 water to measure the chemical potential at 470 K and 0.861 g/cc. I found that the chemical potential was 20%-25% higher than literature values. I believe there may be a problem related to using PME with test-particle insertion.

To test the PME code, I ran a very short TPI run of a water configuration with one molecule deleted. I forced insertions into the hole by using "tpic" mode. I compared the potential energies of these insertions with those obtained by simply taking the difference in potential energy of two md runs, where one contains the missing molecule and one has a hole instead. I ran mdrun in with the "-debug 1" option and added a line of code to tpi.c to print out the insertion locations. I am running the "release-4-5-patches origin/release-4-5-patches branch" of Gromacs 4.5.1.

The LJ and short-range coulomb energies agree well, but the long-range coulomb energy does not. Compare 0.0229272 kJ/mol in TPI mode with -8.771680 kJ/mol in MD mode. Note that the total energy changes are -61.19020 and -69.97481, respectively.

input_files.tar.gz (291 KB) input_files.tar.gz Input files Kevin Daly, 10/08/2010 09:29 AM
tpi_inp_pme.tar.gz (1.41 MB) tpi_inp_pme.tar.gz Input files to test bug fix Kevin Daly, 10/13/2010 06:40 AM
tpi.c (26.9 KB) tpi.c Modified tpi.c that prints out insertion coordinates Kevin Daly, 10/13/2010 04:26 PM

History

#1 Updated by Kevin Daly about 9 years ago

Created an attachment (id=550)
Input files

#2 Updated by Berk Hess about 9 years ago

Could you attach an mdp file for tpi as well?

Berk

#3 Updated by Berk Hess about 9 years ago

I put an optimization in the PME code after the TPI PME implementation
where the grid was not back-fft'ed when forces were not required.
I fixed it:
commit be0a266eb1681c1a3360db32606afdda0174c8eb
I reran and now the energies seem to be correct.

I now trying to find out why tpi with pme does not run in parallel.

Berk

#4 Updated by Berk Hess about 9 years ago

I now also fixed TPI not running in parallel with PME.

Berk

#5 Updated by Kevin Daly about 9 years ago

Thank you for taking the time to fix the LR electrostatics. I ran some tests again, and the LR energies are now of the same order of magnitude, but still differ appreciably (-3.8486 vs. 7.4353 for TPI and brute force, respectively). This difference does not get any smaller if I increase the PME order/grid size or number of wavevectors for Ewald.

I have attached a complete set of input files for the tests as well as a summary of the results. I compiled the release-4.5-patches branch of the Gromacs code with the most recent commit being 104e26658fc6c4dd28867ecb31868b9ea398780c.

#6 Updated by Kevin Daly about 9 years ago

Created an attachment (id=554)
Input files to test bug fix

#7 Updated by Berk Hess about 9 years ago

But you can't directly compare one frame with plain Gromacs.
You seem to be inserting the molecule only once.
The water molecule is placed and oriented randomly in the "cavity",
so you'll never get matching energies.

When fixing the bug I tested the code by on your system by reinserting a water
molecule in the cavity in 10000 times in each of 101 trajectory frames.
The PME energy varies a lot per frame, as I would expect.
The PME energy difference I got is -13.9 with a std.dev. of 4.4.

Note that you should center the water molecule on the oxygen
in the -c file for grompp and the trajectory should have the coordinate
of the oxygen as last coordinate for optimal sampling.

Berk

#8 Updated by Kevin Daly about 9 years ago

(In reply to comment #7)

But you can't directly compare one frame with plain Gromacs.
You seem to be inserting the molecule only once.
The water molecule is placed and oriented randomly in the "cavity",
so you'll never get matching energies.

When fixing the bug I tested the code by on your system by reinserting a water
molecule in the cavity in 10000 times in each of 101 trajectory frames.
The PME energy varies a lot per frame, as I would expect.
The PME energy difference I got is -13.9 with a std.dev. of 4.4.

Note that you should center the water molecule on the oxygen
in the -c file for grompp and the trajectory should have the coordinate
of the oxygen as last coordinate for optimal sampling.

Berk

In the particular tests that I ran I actually modified tpi.c to print out the randomly generated coordinates of all four particles in the inserted molecule. With ld_seed = 1034, I found them to be:
0.88429457 4.09253597 1.03514338
0.85779530 4.17442083 0.99384421
0.97501451 4.07984591 1.00865090
0.89276773 4.10166550 1.02619910
for oxygen, hydrogen1, hydrogen2, and the virtual site respectively. I then took these same coordinates and placed them in both the configuration and trajectory files of the N+1 system. Since I used the exact coordinates of the insertion, I would expect to get good agreement, and the LJ, LJ_corr, and Coul_SR energies do agree quite well. It is only the Coul_LR energy that seems to be different.

I have attached the modified version of tpi.c.

-Kevin

#9 Updated by Kevin Daly about 9 years ago

Created an attachment (id=555)
Modified tpi.c that prints out insertion coordinates

#10 Updated by Berk Hess about 9 years ago

Ah, after a long search I found out the charges of the inserted molecule
were put on the grid. This was due to a change in the pme code just before
the 4.5 release. I fixed it along with a periodicity issue in the TPI+PME
combination:
commit 43ffe2f7c5171b9970bb179f17b2570d78a2232d

Berk

Also available in: Atom PDF