Re: Kb, Ktheta values for TIP3P water model

From: Maxim Belkin (mbelkin_at_ks.uiuc.edu)
Date: Tue Aug 12 2014 - 23:11:43 CDT

Viswanath,

In NAMD, 1 rigid water molecule adds "2" to "Info: * BONDS", and "3" to "Info: * RIGID BONDS". So, NAMD simulations of 200 rigid water molecules should report 200 angles (Info: * Angles), 400 bonds (Info: * Bonds), and 600 rigid bonds (Info: * rigid bonds). In this respect, seeing 600 rigid bonds for 200 water molecules in NAMD log is fine.

Did you notice that your system is not neutral? This is not great.

> TOTAL CHARGE = 18.0571 e

Perhaps, it would be better to rebuild your system. You can explicitly remove (or comment out) hydrogen-hydrogen bond from your topology file.

Another question: can you find the following lines in your log file?

Info: RIGID BONDS TO HYDROGEN : ALL
Info: ERROR TOLERANCE : 1e-08
Info: MAX ITERATIONS : 100
Info: RIGID WATER USING SETTLE ALGORITHM

Maxim

On Aug 12, 2014, at 12:19 PM, Viswanath Pasumarthi <v.pasumarthi_at_iitg.ernet.in> wrote:

> Yes Kenno, there is a warning message which I have failed to notice when
> force constant values for O-H bond and H-O-H angle are taken to be zero.
>
> Warning: Ignored 400 bonds with zero force constants.
> Warning: Will get H-H distance in rigid H2O from H-O-H angle.
> Warning: Ignored 200 angles with zero force constants.
>
> Moving on, as Maxim suggested, on checking the variation of length of O-H
> bond with force constants, kb = 1,2,5,10,100,450, the bond lengths varied
> from 5% to 15% on average for 10 randomly selected bonds. This only
> reinforces the failure of rigidBonds all constraint. But the second
> warning message posted above indicates that NAMD identified water as a
> rigid molecule.
>
> The following is an excerpt from output log file for a system involving
> ionic liquid+water mixture. 20 IL molecules + 200 H2O molecules were
> considered. It can be seen that 900 rigid bonds are identified. Each IL
> molecule contains 15 covalent bonds involving hydrogen, which makes 300
> rigid bonds corresponding to IL. The rest 600, contributes from 200 water
> molecules, which implies even H-H is considered as an explicit bond for
> water molecule. Is it alright and does it mean the bond parameters for H-H
> should also be provided?
>
> Info: STRUCTURE SUMMARY:
> Info: 1400 ATOMS
> Info: 1180 BONDS
> Info: 1600 ANGLES
> Info: 1660 DIHEDRALS
> Info: 0 IMPROPERS
> Info: 0 CROSSTERMS
> Info: 0 EXCLUSIONS
> Info: 900 RIGID BONDS
> Info: 3297 DEGREES OF FREEDOM
> Info: 700 HYDROGEN GROUPS
> Info: 4 ATOMS IN LARGEST HYDROGEN GROUP
> Info: 700 MIGRATION GROUPS
> Info: 4 ATOMS IN LARGEST MIGRATION GROUP
> Info: TOTAL MASS = 11990.3 amu
> Info: TOTAL CHARGE = 18.0571 e
> Info: MASS DENSITY = 0.73744 g/cm^3
> Info: ATOM DENSITY = 0.0518519 atoms/A^3
>
>
> PSF:
> 801 MOL 1 UNK O OW -0.928231 15.9994 0
> 802 MOL 1 UNK H HT 0.464116 1.00794 0
> 803 MOL 1 UNK H HT 0.464116 1.00794 0
>
>
> BONDS
> !atom type Kb b0
> !TIP3P WATER
> OW HT 450.0000 0.9572 !
>
> ANGLES
> !atom types Ktheta Theta0 Kub S0
> !TIP3P Water
> HT OW HT 55.0000 104.52 !
>
> NONBONDED
> !atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
> !TIP3P Water
> HT 0.000000 -0.046000 0.224500 !
> OW 0.000000 -0.152100 1.768200 !
>
>
>
>> On 08/11/2014 01:56 PM, Viswanath Pasumarthi wrote:
>>> considering a value of k=0 mean absence of spring (physical
>>> bond) between the respective atoms. It would have been easier to
>>> understand this mistake had there been an error message pointing to the
>>> choice of k value as zero.
>>
>> But there *is* a message, and if you would have noticed/posted it, we
>> wouldn't have needed to have this long, marginally constructive
>> discussion. Look in your NAMD output for the following:
>>
>> Warning: Ignored .* bonds with zero force constants.
>> Warning: Ignored .* angles with zero force constants.
>>
>> Where .* is an integer.
>>
>>> However, using rigidBonds constraints still
>>> required bond constraint values to be specified in the parameter file.
>>> This, as Kenno cautioned has turned into a flexible simulation, as
>>> covalent bonds involving hydrogen became more flexible at lower k values
>>> (1,2,5,10) and functioning decent at the values reported in modified
>>> TIP3P
>>> model.
>>
>> I have to repeat that I tried and failed to reproduce this; my rigidBonds
>> simulations behave correctly for all force constants (even absurd values
>> such as 0.0001 or -10000) except zero. Most probably, there is an issue
>> with your psf and/or rigidBonds is not applied correctly. Which I have
>> been saying from the very beginning!
>>
>> You probably need to try to use more detailed language (for example
>> "distances" or "force constants" rather than "values") when describing
>> your problem, and attach relevant input and output files to your posts,
>> else we'll never get anywhere.
>>
>>

This archive was generated by hypermail 2.1.6 : Thu Dec 31 2015 - 23:21:06 CST