Project

General

Profile

Bug #1115

Using PME-Switch with Thermodynamic Integration

Added by Thomas Purcell almost 5 years ago. Updated over 3 years ago.

Status:
Closed
Priority:
Normal
Category:
mdrun
Target version:
Affected version - extra info:
4.6, 4.5
Affected version:
Difficulty:
uncategorized
Close

Description

I am using Gromacs to run a TI calculation to study an amino acid mutation. I found that the Coulomb-short range energy (Coulomb-SR) term when TI was turned on with a lambda value of 0.0 was not the same as when TI was turned off. In order to determine if it was a topology problem I ran the same system with a standard topology file of state A with the same coordinate file with the dummy atoms removed. The Coulomb-SR interaction for the standard topology agreed with the results from the dual topology when TI was turned off. After various tests with the TI parameters it was discovered the problem was with the PME-Switch method. When PME or Cut-off was used the Coulomb-short range energy agreed.

Attached is a .tar.gz with files to show the energy difference between PME and PME-switch, and force field files that have been modified to handle for the new interactions for the dummy atoms.

TI_PME-Switch_problem.tar.gz (552 KB) TI_PME-Switch_problem.tar.gz Directory of files comparing the standard conformation with the TI one using PME and PME-Switch Thomas Purcell, 01/11/2013 04:38 PM
charmm27_TI_modified.tar.gz (63.5 KB) charmm27_TI_modified.tar.gz Modfied Charmm forcefield files (adds the new atom types and interactions) Thomas Purcell, 01/11/2013 04:38 PM
0001-Fixes-the-Function-type-U-B-not-implemented-in-ip_pe.patch (1.13 KB) 0001-Fixes-the-Function-type-U-B-not-implemented-in-ip_pe.patch Magnus Lundborg, 03/15/2013 04:00 PM

Associated revisions

Revision 780247b0 (diff)
Added by Magnus Lundborg over 4 years ago

Fixes "Function type U-B not implemented in ip_pert."

Fixes this bug, occuring when running free energy calculations using
the CHARMM ff. Refs #1115.

Change-Id: Ia3becd61eca93447c936bdb3297570684ae29a28

Revision 2ed41ed5 (diff)
Added by Magnus Lundborg over 4 years ago

Fixes "Function type U-B not implemented in ip_pert."

Fixes this bug, occuring when running free energy calculations using
the CHARMM ff. Refs #1115.

Change-Id: Ia3becd61eca93447c936bdb3297570684ae29a28

History

#1 Updated by Michael Shirts almost 5 years ago

  • Category set to mdrun
  • Assignee set to Michael Shirts

#2 Updated by Michael Shirts almost 5 years ago

Not sure if the error is in 4.6 yet because I get:

"Function type U-B not implemented in ip_pert".

Will have to try a workaround.

#3 Updated by Michael Shirts almost 5 years ago

Error does appear to exist with this topology in 4.6 (I just commented out all the [ angles ] and got around the previous error, though it should be looked at as well -- why is it not defined for ip_pert in 4.6?)

Will investigate, though won't be able to until tomorrow. Difference is about 0.2 kJ/mol, i.e. -635.441643 vs -635.659431, enough to pay attention to.

#4 Updated by Michael Shirts almost 5 years ago

  • Target version set to 4.6
  • Affected version - extra info set to 4.6, 4.5

#5 Updated by Thomas Purcell almost 5 years ago

Hi,

For the i_pert problem here is the work around I use

insert into line 111 of topsort.c
case F_CMAP:
bPert = FALSE;
gmx_warning("FEP or TI with charmm27 work only if the CMAP term is unchanged between A and B.");

That should work for this system, I am not sure if it does overall.

I also ran into a problem with the UB terms not interpolating if this is still a problem in version 6 I can attach a file with the code to fix that problem

#6 Updated by Erik Lindahl almost 5 years ago

  • Target version changed from 4.6 to 4.6.1

This was unfortunately reported long after our deadlines for 4.6, so it won't be a blocker for that release. Bumping to 4.6.1.

#7 Updated by Michael Shirts almost 5 years ago

Agreed on the timing.

#8 Updated by Michael Shirts almost 5 years ago

OK, got a chance to look at this. Looking at the equil.mdp file, I see

coulombtype = PME-switch
rlist = 1.4
rcoulomb = 0.9
vdwtype = switch

When I look at the mdp.out, I see:

coulombtype = PME-switch
rcoulomb-switch = 0
rcoulomb = 0.9

So the rcoulomb-switch is being set to zero, leading to a very odd potential. I wonder if this is the source of the error?

When I set rcoulomb switch to something reasonable, I get good agreement between the energies from PME-switch with free energies on and off.

Thomas, when you make this change, do you also get correct results again?

#9 Updated by Thomas Purcell almost 5 years ago

I do get better agreement, but they should be an exact match because the interactions they measure are the same if init_lambda = 0, and TI is turned off.

With a more reasonable rcoulomb-switch value I get an agreement to the third decimal place.

#10 Updated by Michael Shirts almost 5 years ago

Thomas Purcell wrote:

I do get better agreement, but they should be an exact match because the interactions they measure are the same if init_lambda = 0, and TI is turned off.

With a more reasonable rcoulomb-switch value I get an agreement to the third decimal place.

OK, this sounds like what I was seeing -- almost complete matching, but not quite exact.

However, I appear to be seeing that the results match up through five digits in 4.6. If this is the case, then it will be much easier to patch, since the solution exists already, and it just needs to be backported. If it's not to hard to check if the release 4.6 code solves this problem (it does for me, but independent verification is useful) then I can go through and identify the code that fixed it.

#11 Updated by Michael Shirts almost 5 years ago

Also, what value of rcoulomb-switch are you using? In most cases, a difference of 0.01 or 0.02 between rcoulomb and rcoulomb-switch is an appropriate amount. The switch should be quite short.

I agree that the values should eventually be exact.

Michael Shirts wrote:

Thomas Purcell wrote:

I do get better agreement, but they should be an exact match because the interactions they measure are the same if init_lambda = 0, and TI is turned off.

With a more reasonable rcoulomb-switch value I get an agreement to the third decimal place.

OK, this sounds like what I was seeing -- almost complete matching, but not quite exact.

However, I appear to be seeing that the results match up through five digits in 4.6. If this is the case, then it will be much easier to patch, since the solution exists already, and it just needs to be backported. If it's not to hard to check if the release 4.6 code solves this problem (it does for me, but independent verification is useful) then I can go through and identify the code that fixed it.

#12 Updated by Thomas Purcell almost 5 years ago

I am using an rcoulomb-switch = 0.89 in version 4.5.5.

I got a segmentation fault when I tried using version 4.6-beta3. I am going to try to fix the error and see what I get.

Michael Shirts wrote:

Also, what value of rcoulomb-switch are you using? In most cases, a difference of 0.01 or 0.02 between rcoulomb and rcoulomb-switch is an appropriate amount. The switch should be quite short.

I agree that the values should eventually be exact.

Michael Shirts wrote:

Thomas Purcell wrote:

I do get better agreement, but they should be an exact match because the interactions they measure are the same if init_lambda = 0, and TI is turned off.

With a more reasonable rcoulomb-switch value I get an agreement to the third decimal place.

OK, this sounds like what I was seeing -- almost complete matching, but not quite exact.

However, I appear to be seeing that the results match up through five digits in 4.6. If this is the case, then it will be much easier to patch, since the solution exists already, and it just needs to be backported. If it's not to hard to check if the release 4.6 code solves this problem (it does for me, but independent verification is useful) then I can go through and identify the code that fixed it.

#13 Updated by Michael Shirts almost 5 years ago

Try the full release, not 4.6-beta3.

Thomas Purcell wrote:

I am using an rcoulomb-switch = 0.89 in version 4.5.5.

I got a segmentation fault when I tried using version 4.6-beta3. I am going to try to fix the error and see what I get.

Michael Shirts wrote:

Also, what value of rcoulomb-switch are you using? In most cases, a difference of 0.01 or 0.02 between rcoulomb and rcoulomb-switch is an appropriate amount. The switch should be quite short.

I agree that the values should eventually be exact.

Michael Shirts wrote:

Thomas Purcell wrote:

I do get better agreement, but they should be an exact match because the interactions they measure are the same if init_lambda = 0, and TI is turned off.

With a more reasonable rcoulomb-switch value I get an agreement to the third decimal place.

OK, this sounds like what I was seeing -- almost complete matching, but not quite exact.

However, I appear to be seeing that the results match up through five digits in 4.6. If this is the case, then it will be much easier to patch, since the solution exists already, and it just needs to be backported. If it's not to hard to check if the release 4.6 code solves this problem (it does for me, but independent verification is useful) then I can go through and identify the code that fixed it.

#14 Updated by Thomas Purcell almost 5 years ago

Still getting segmentation fault

Michael Shirts wrote:

Try the full release, not 4.6-beta3.

Thomas Purcell wrote:

I am using an rcoulomb-switch = 0.89 in version 4.5.5.

I got a segmentation fault when I tried using version 4.6-beta3. I am going to try to fix the error and see what I get.

Michael Shirts wrote:

Also, what value of rcoulomb-switch are you using? In most cases, a difference of 0.01 or 0.02 between rcoulomb and rcoulomb-switch is an appropriate amount. The switch should be quite short.

I agree that the values should eventually be exact.

Michael Shirts wrote:

Thomas Purcell wrote:

I do get better agreement, but they should be an exact match because the interactions they measure are the same if init_lambda = 0, and TI is turned off.

With a more reasonable rcoulomb-switch value I get an agreement to the third decimal place.

OK, this sounds like what I was seeing -- almost complete matching, but not quite exact.

However, I appear to be seeing that the results match up through five digits in 4.6. If this is the case, then it will be much easier to patch, since the solution exists already, and it just needs to be backported. If it's not to hard to check if the release 4.6 code solves this problem (it does for me, but independent verification is useful) then I can go through and identify the code that fixed it.

#15 Updated by Michael Shirts almost 5 years ago

Thomas Purcell wrote:

Still getting segmentation fault

Strange. Can you run ANYTHING with 4.6, or just not these files? I'm not having any issues running. Perhaps it's an installation problem?

#16 Updated by Thomas Purcell almost 5 years ago

It's probably an installation error. I can look at it tomorrow when I have more time.

#17 Updated by Thomas Purcell almost 5 years ago

I am not entirely sure what was causing the segmentation faults with mdrun on the computer I was using earlier I am still getting them so I shall do a more thorough check later. Once I installed it on a different computer I got it to work by commenting out all the angles like you did. The results agreed to 5 significant figures like you had.

#18 Updated by Magnus Lundborg almost 5 years ago

This patch fixes the "Function type U-B not implemented in ip_pert" problem (for me). Please test it.

#19 Updated by Thomas Purcell over 4 years ago

Sorry I was on vacation for the past 2 weeks, I shall test it by the end of the week.

#20 Updated by Thomas Purcell over 4 years ago

The fix works for me too.

#21 Updated by Mark Abraham over 4 years ago

  • Status changed from New to Accepted
  • Assignee changed from Michael Shirts to Magnus Lundborg
  • Target version deleted (4.6.1)
  • Affected version set to 4.6

Magnus, can you upload your fix to gerrit, please?

#22 Updated by Magnus Lundborg over 4 years ago

I have now uploaded my fix (https://gerrit.gromacs.org/#/c/2346/). I still don't think the whole issue is solved as I just tried to fix a small part of the problem. But if Thomas's comment meant that it fixed the bug for him this can be closed.

#23 Updated by Thomas Purcell over 4 years ago

Magnus Lundborg wrote:

I have now uploaded my fix (https://gerrit.gromacs.org/#/c/2346/). I still don't think the whole issue is solved as I just tried to fix a small part of the problem. But if Thomas's comment meant that it fixed the bug for him this can be closed.

It fixed the issue with the Function type U-B not implemented in ip_pert

I still did not get exactly the same values using PME-Switch

#24 Updated by Mark Abraham over 4 years ago

  • Status changed from Accepted to In Progress

https://gerrit.gromacs.org/#/c/2340/ co-incidentally fixes the aspect of the problem where rcoulomb-switch defaults to a silly value of 0 which is silently accepted.

#25 Updated by Mark Abraham over 4 years ago

Unfortunately I don't know enough about how PME-switch or TI should work to be of further help.

#26 Updated by Rossen Apostolov over 3 years ago

  • Assignee changed from Magnus Lundborg to Michael Shirts
  • Target version set to 4.6.x

Michael, Thomas: what is the current situation with the energy differences issue?

#27 Updated by Erik Lindahl over 3 years ago

  • Status changed from In Progress to Resolved

This has been fully fixed by the recent free energy and modifier patches that went into 4.6 - I get perfectly matching results for some of the attached files.

However, if there are still remaining issues, please open a new issue (mark this one as related), and make sure the bug report is as simple and specific as possible (only one uploaded system, possible with two alternative mdp files).

#28 Updated by Erik Lindahl over 3 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF