Project

General

Profile

Bug #284

couple-intramol not detected in .tpr files by gmxcheck or gmxdump

Added by Chris Neale almost 11 years ago. Updated almost 11 years ago.

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

Description

When free_energy=yes and two .tpr files are created with the following difference:

couple-intramol = yes vs. couple-intramol = no

This couple-intramol parameter difference is not picked up by gmxcheck -s1 -s2

Further, gmxdump -s outputs a list of .mdp parameters that does not include couple-intramol:

...
sa_surface_tension = 2.092
DispCorr = EnerPres
free_energy = yes
init_lambda = 1
sc_alpha = 0
sc_power = 1
sc_sigma = 0.3
delta_lambda = 0
nwall = 0
...

At this point it appears to be a trivial case, as long as this is simply an output issue and not something that doesn't make it into the .tpr file.

Thanks,
Chris.

History

#1 Updated by Berk Hess almost 11 years ago

couple-intramol is a grompp processing parameter
that functions similar to the constraint option.
It modifies your bonded and non-bonded interactions
and mdrun does not need to know about couple-intramol,
so it is not stored in the tpr file.

This might be slightly inconvenient for checking
your tpr files, but it is not a bug.

You will see a lot of differences in interactions
between couple-intramol on and off.

Berk

#2 Updated by Chris Neale almost 11 years ago

Thanks Berk,

The fact that the flag itself doesn't make it into the .tpr file doesn't concern me too much. I am more concerned about the fact that in fact I do not see any differences in interactions between couple-intramol on and off, even for lambda=1.0

That the couple-intramol flag does not make it into the .tpr file is not even inconvenient. That the differing interactions are not picked up by

While using:
; Free energy control stuff
free_energy = yes
init_lambda = 1.00
delta_lambda = 0
sc_alpha = 0.0
sc-power = 1.0
sc-sigma = 0.3
couple-moltype =
couple-lambda0 = vdw-q
couple-lambda1 = vdw-q

And:
A) couple-intramol = no; to generate nocoupleintramol.tpr
B) couple-intramol = yes; to generate yescoupleintramol.tpr

$ gmxcheck -s1 nocoupleintramol.tpr -s2 yescoupleintramol.tpr
...
Reading file nocoupleintramol.tpr, VERSION 4.0.3 (double precision)
Reading file yescoupleintramol.tpr, VERSION 4.0.3 (double precision)
comparing inputrec
comparing top
comparing idef
comparing ilist BONDS
comparing ilist G96BONDS
comparing ilist MORSE
comparing ilist CUBICBONDS
comparing ilist CONNBONDS
comparing ilist HARMONIC
comparing ilist FENEBONDS
comparing ilist TABBONDS
comparing ilist TABBONDSNC
comparing ilist ANGLES
comparing ilist G96ANGLES
comparing ilist CROSS_BOND_BOND
comparing ilist CROSS_BOND_ANGLE
comparing ilist UREY_BRADLEY
comparing ilist QANGLES
comparing ilist TABANGLES
comparing ilist PDIHS
comparing ilist RBDIHS
comparing ilist FOURDIHS
comparing ilist IDIHS
comparing ilist PIDIHS
comparing ilist TABDIHS
comparing ilist LJ14
comparing ilist COUL14
comparing ilist LJC14_Q
comparing ilist LJC_NB
comparing ilist LJ_SR
comparing ilist BHAM
comparing ilist LJ_LR
comparing ilist BHAM_LR
comparing ilist DISPCORR
comparing ilist COUL_SR
comparing ilist COUL_LR
comparing ilist RF_EXCL
comparing ilist COUL_RECIP
comparing ilist DPD
comparing ilist POLARIZATION
comparing ilist WATERPOL
comparing ilist THOLE
comparing ilist POSRES
comparing ilist DISRES
comparing ilist DRVIOL
comparing ilist ORIRES
comparing ilist ORDEV
comparing ilist ANGRES
comparing ilist ANGRESZ
comparing ilist DIHRES
comparing ilist DIHVIOL
comparing ilist CONSTR
comparing ilist CONSTRNC
comparing ilist SETTLE
comparing ilist VSITE2
comparing ilist VSITE3
comparing ilist VSITE3FD
comparing ilist VSITE3FAD
comparing ilist VSITE3OUT
comparing ilist VSITE4FD
comparing ilist VSITE4FDN
comparing ilist VSITEN
comparing ilist COM_PULL
comparing ilist EQM
comparing ilist EPOT
comparing ilist EKIN
comparing ilist ETOT
comparing ilist ECONS
comparing ilist TEMP
comparing ilist PRES
comparing ilist DV/DL
comparing ilist DK/DL
comparing ilist DG/DL_CON
comparing atoms
comparing block cgs
comparing block mols
comparing blocka excls
comparing flags
comparing box
comparing box_rel
comparing boxv
comparing x
comparing v

gcq#8: "Don't Push Me, Cause I'm Close to the Edge" (Tricky)

#########

$ gmxdump -s nocoupleintramol.tpr > no.dump
$ gmxdump -s yescoupleintramol.tpr > yes.dump

$ diff no.dump yes.dump
1c1
< nocoupleintramol.tpr:
---

yescoupleintramol.tpr:

#########

And then also note that this segment of the gmxdump output appears to confirm to me that there is indeed a difference between the A and B state charge value that makes it into the .tpr

...
moltype (0):
name="DPN"
atoms:
atom (23):
atom[ 0]={type= 0, typeB= 0, ptype= Atom, m= 1.50350e+01, q= 4.00000e-01, mB= 1.50350e+01, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom[ 1]={type= 0, typeB= 0, ptype= Atom, m= 1.50350e+01, q= 4.00000e-01, mB= 1.50350e+01, qB= 0.00000e+00, resnr= 0, atomnumber= 1}
atom[ 2]={type= 0, typeB= 0, ptype= Atom, m= 1.50350e+01, q= 4.00000e-01, mB=
...

where this difference does extend to atoms that are separated by more than 3 bonds.

##########

However, I do see that the lambda term makes it into the .tpr file based on this output from gmxdump:

...
header:
bIr = present
bBox = present
bTop = present
bX = present
bV = present
bF = not present
natoms = 29870
step = 0
t = 0.000000e+00
lambda = 1.000000e+00
...

##########

And here is the most worrying thing for me:

$diff nocoupleintramol.tpr yescoupleintramol.tpr
(it returns no difference)

while 'diff' picks up a difference when I set free_energy=no:
$ diff nofreeenergy.tpr yescoupleintramol.tpr
Files nofreeenergy.tpr and yescoupleintramol.tpr differ

##########

Given that the .tpr files are identical (not withstanding the fact that the mdout.mdp files generated by grompp differ), it is perhaps not surprising that the energies from mdrun are also identical from couple-intramol=yes or =no.

##########

Perhaps I am misunderstanding something about

couple-moltype =
couple-lambda0 = vdw-q
couple-lambda1 = vdw-q

which I actually did not include in my .mdp files, but copied from the generated mdout.mdp for the same of completeness.

From the manual it appears to me that these options might be used instead of directly modifying a topology and adding a B-state that has, for example, no charges.

Thanks again,
Chris.

#3 Updated by Chris Neale almost 11 years ago

Here are the energy values to which I was referring:

couple-intramol=no; init_lambda=1.0
Energies (kJ/mol)
Angle Proper Dih. Ryckaert-Bell. LJ-14 Coulomb-14
2.91166e+03 5.90834e+02 1.07820e+03 7.70226e+02 6.09619e+03
LJ (SR) LJ (LR) Disper. corr. Coulomb (SR) Coul. recip.
4.94951e+04 -1.84944e+03 -6.38472e+02 -3.22179e+05 -7.03149e+04
Potential Kinetic En. Total Energy Temperature Pressure (bar)
-3.34039e+05 5.67730e+04 -2.77266e+05 3.00622e+02 6.78592e+01
dVpot/dlambda dEkin/dlambda dG/dl constr. Cons. rmsd () Cons.2 rmsd ()
1.29828e+03 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

couple-intramol=yes; init_lambda=1.0
Energies (kJ/mol)
Angle Proper Dih. Ryckaert-Bell. LJ-14 Coulomb-14
2.91166e+03 5.90834e+02 1.07820e+03 7.70226e+02 6.09619e+03
LJ (SR) LJ (LR) Disper. corr. Coulomb (SR) Coul. recip.
4.94951e+04 -1.84944e+03 -6.38472e+02 -3.22179e+05 -7.03149e+04
Potential Kinetic En. Total Energy Temperature Pressure (bar)
-3.34039e+05 5.67730e+04 -2.77266e+05 3.00622e+02 6.78654e+01
dVpot/dlambda dEkin/dlambda dG/dl constr. Cons. rmsd () Cons.2 rmsd ()
1.29828e+03 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

###

couple-intramol=no; init_lambda=0.5
Energies (kJ/mol)
Angle Proper Dih. Ryckaert-Bell. LJ-14 Coulomb-14
2.91166e+03 5.90834e+02 1.07820e+03 7.70226e+02 6.15951e+03
LJ (SR) LJ (LR) Disper. corr. Coulomb (SR) Coul. recip.
4.94951e+04 -1.84944e+03 -6.38472e+02 -3.22416e+05 -7.07897e+04
Potential Kinetic En. Total Energy Temperature Pressure (bar)
-3.34688e+05 5.67662e+04 -2.77922e+05 3.00586e+02 3.49030e+01
dVpot/dlambda dEkin/dlambda dG/dl constr. Cons. rmsd () Cons.2 rmsd ()
1.29828e+03 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

couple-intramol=yes; init_lambda=0.5
Energies (kJ/mol)
Angle Proper Dih. Ryckaert-Bell. LJ-14 Coulomb-14
2.91166e+03 5.90834e+02 1.07820e+03 7.70226e+02 6.15951e+03
LJ (SR) LJ (LR) Disper. corr. Coulomb (SR) Coul. recip.
4.94951e+04 -1.84944e+03 -6.38472e+02 -3.22416e+05 -7.07897e+04
Potential Kinetic En. Total Energy Temperature Pressure (bar)
-3.34688e+05 5.67662e+04 -2.77922e+05 3.00586e+02 3.49066e+01
dVpot/dlambda dEkin/dlambda dG/dl constr. Cons. rmsd () Cons.2 rmsd ()
1.29828e+03 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00

#4 Updated by Chris Neale almost 11 years ago

Note that when I create a DPN molecule that is identical to DPC (no defined B-state) and utilize "couple-moltype = DPN", I do see the type of difference from "gmxcheck -s1 nocoupleintramol.tpr -s2 yescoupleintramol.tpr" that I expect.

Could it be that couple-intramol is nonfunctional when the user explicitly defines the B-state in the .itp file and then does not utilize couple-moltype?

Thanks,
Chris.

#5 Updated by Chris Neale almost 11 years ago

Oops, I see that this is actually explained in the manual, which I did read but obviously not closely enough.

"All intra-molecular non-bonded interactions for moleculetype couple-moltype are replaced by exclusions and explicit pair interactions. In this manner the decoupled state of the molecule corresponds to the proper vacuum state without periodicity effects."

You can close this ticket again.

Sorry for wasting your time.
Chris.

#6 Updated by Chris Neale almost 11 years ago

I suggest that grompp throw a fatal error upon the explicit declaration of couple-intramol when couple-moltype is undefined.

I realize that this may not be possible, but if it is, it would help users avoid thinking that they are getting couple-intramol=no when really they are not since they specified the B-state explicitly in the itp file.

Chris.

#7 Updated by Berk Hess almost 11 years ago

I agree that it would be nice if grompp could give a fatal error.
But on the other hand I like the feature for these mdp parameters
"blocks" that if you a block off, couple-moltype=no,
all parameters are ignored and you do not get a fatal error.

Also the mdp manual clearly states that couple-intramol
is linked to couple-moltype.

Berk

Also available in: Atom PDF