Project

General

Profile

Bug #1318

Diverging behavior of endstate with and without FE code.

Added by Manuel Melo over 4 years ago. Updated about 4 years ago.

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

Description

Hi,

I noticed in 4.6.3 an odd behavior, on an NPT simulation of pre-equilibrated boxes of different solvents:

If I turn on the free-energy code (free-energy=yes) and set the system to a fully coupled endstate (init-lambda=1; couple-moltype=system; couple-lambda1=vdw-q; couple-lambda0=vdw) the resulting run diverges significantly from a simulation run with free-energy=no. Namely, the pressure coupling just causes the system to continuously expand.

The exact same mdps run fine in 4.5.7, with and without the free-energy code. Below is the mdp I'm using (including the FE section; for plain runs I simply omitted that section).

For reference, I have observed this behavior with spc water and gromos53a6 propanol (other 53a6 solvents also seemed to be affected but I did not look into that thoroughly). The mdp and attached tprs and plot are for propanol.
The bug also manifests itself when using A-B topologies in the itp rather than in the mdp. Version 4.6.2 also seemed to be affected.

Please let me know what other feedback I can provide.

Cheers,
Manel

#####################################
integrator = md
tinit = 0.0
dt = 0.002
nsteps = 1500000
comm-mode = Linear
nstcomm = 100
comm-grps = System

nstxout = 1000
nstvout = 0
nstfout = 1000
nstlog = 5000
nstenergy = 5000
nstxtcout = 0

nstlist = 10
ns_type = grid
pbc = xyz
rlist = 0.9
rlistlong = 1.4

coulombtype = reaction-field
rcoulomb = 1.4
epsilon_r = 1
epsilon_rf = 19.4
vdw-type = Cut-off
rvdw = 1.4

Tcoupl = V-rescale
tc-grps = system
tau_t = 0.1
ref_t = 298
nsttcouple = 10
Pcoupl = berendsen
Pcoupl_type = isotropic
tau_p = 1
compressibility = 1e-5
ref_p = 1.0
nstpcouple = 10

gen_vel = yes
gen_temp = 298
gen_seed = -1

constraints = hbonds

;;FREE-ENERGY;;;;;;;;;;;;;;;;;;;;;;;;
free-energy = yes
init-lambda = 1
couple-moltype = system
couple-lambda1 = vdw-q
couple-lambda0 = vdw
nstdhdl = 0

#####################################

FEbug.png (8.43 KB) FEbug.png Manuel Melo, 08/09/2013 12:25 PM
tprs.tgz (382 KB) tprs.tgz Manuel Melo, 08/09/2013 12:25 PM
FE-propanol.tgz (553 KB) FE-propanol.tgz Manuel Melo, 09/28/2013 11:12 PM

Associated revisions

Revision 6730ca88 (diff)
Added by Berk Hess about 4 years ago

Fixes reaction field free energy bug

Version 4.6 introduced an error in the reaction-field correction
force term for perturbed interactions. A factor r_softcore^2 was
missing. The force calculation code is now slightly reorganized
and comments added to avoid such issues in the future.
Fixes #1318

Change-Id: I9105139f8975495c323008ce202cde517a69281a

History

#1 Updated by Manuel Melo over 4 years ago

Possibly related to bug #1315.

#2 Updated by Michael Shirts about 4 years ago

  • Category set to mdrun
  • Assignee set to Michael Shirts
  • Target version set to 4.6.4

#3 Updated by Michael Shirts about 4 years ago

Can you add .mpd/.top/.gro files so that I can try to tweak some things as well as understand more completely what is going on?

Thanks!

#4 Updated by Manuel Melo about 4 years ago

Here you go. The directory tree in the .tgz should be self-evident.

Cheers, and thanks for looking into this!

#5 Updated by Michael Shirts about 4 years ago

I've traced it down to the force (but not the energy) in the inner loops being calculated differently between the two simulations, but slow going from here.

#6 Updated by Michael Shirts about 4 years ago

It appears to be a problem with the free energy kernel running reaction field, not twin range or any other part.

I get agreement when switch to

rlist = 1.0

coulombtype = PME
rcoulomb = 1.0
vdw-type = Switch
rvdw = 0.9
rvdw-switch = 0.88

But when I switch coulombtype to
coulombtype = reaction-field

Then I get disagreement between free energy and non-free energy.

Will see how much further I can get tonight.

#7 Updated by Michael Shirts about 4 years ago

  • Status changed from New to Fix uploaded

OK, I believe I found the error. When the free energy inner loops were written for 4.6, there was an algebra mistake. A checked in fix should resolve the problem.

#8 Updated by Manuel Melo about 4 years ago

Cheers, and thanks for your time.

#9 Updated by Michael Shirts about 4 years ago

Manuel Melo wrote:

Cheers, and thanks for your time.

Thanks for providing useful and important debugging information.

#10 Updated by Berk Hess about 4 years ago

  • Status changed from Fix uploaded to Resolved
  • % Done changed from 0 to 100

#11 Updated by Mark Abraham about 4 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF