[gmx-developers] Re:Re:Re:Re:Re:Re: Energy Minimization

[Leontyev Igor]

> Igor Leontyev wrote:
>>> I just noticed that settle for the coordinates (the forces were correct)
>>> does not support partially frozen water molecules.
>>> Try replacing in src/mdlib/mdatom.c:
>>>      if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] &&
>>> opts->nFreeze[g][ZZ])
>>>    /* Set the mass of completely frozen particles to ALMOST_ZERO iso 0
>>>     * to avoid div by zero in lincs or shake.
>>>     * Note that constraints can still move a partially frozen particle.
>>>     */
>>>    md->invmass[i]    = ALMOST_ZERO;
>>>      else
>>> by
>>>      if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] &&
>>> opts->nFreeze[g][ZZ]) {
>>>    /* Set the mass of completely frozen particles to ALMOST_ZERO iso 0
>>>     * to avoid div by zero in lincs or shake.
>>>     * Note that constraints can still move a partially frozen particle.
>>>     */
>>>    md->massT[i] = 1/ALMOST_ZERO;
>>>    md->invmass[i]    = ALMOST_ZERO;
>>>      } else
>>>
>>> Berk
>>>
>> Minimization with the code replacement produced exactly! the same output.
>> Igor Leontyev

> Ah, there might be another problem as well.
> You have some completely frozen waters and some partially frozen waters,
> right?
> The current SETTLE code assumes the masses of all waters are identical
> and does not check for this.
>
> If you use constraints=all-bonds, you can use the mdp option:
> define = -DFLEXIBLE
> which will make the water model flexible, but constraints=all-bonds will
> replace
> the bonds by constraints, which will then use LINCS iso SETTLE.
> Berk
>

I use to try the trick. However,  "constraints=all-bonds" with flexible 
water
result on no optimization at all. I don't know why in my case an orientation 
of optimizing water molecules is not changing. AmberFF tip3p water with
added ability of flexibility is used:
---------------------------------
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.09572 345000 0.09572 345000
1 3 1 0.09572 345000 0.09572 345000
[ angles ]
; i j k funct angle force.c.
2 1 3 1 105.4 383 105.4 383
#else
[ settles ]
; OW funct doh dhh
1 1 0.09572 0.15139
[ exclusions ]
1 2 3
2 1 3
3 1 2
#endif
---------------------------------

Use of "constraints=hbonds" (instead of "all-bonds") with the flexible water 
does the orientational optimization much better. However, angles of water 
molecules are distorted in the final structure.
Igor