[gmx-users] Can drude shell positions be minimized with keeping atoms positions restrained?

[Jones De Andrade]

Hi everybody.

Since my first email seemed to be not following the list's nettiquete 
rules (sorry for that), I'm trying again.

First, I was not successful in solving this issue on gromacs versions 
2018.4, 2019.6 and 2020.1: so, it is consistent among them.

What we are trying to do is to make the energy optimization only for 
the Drude shells of the molecules we are using. At this point we tried 
both steepest-descent and conjugated-gradient, but when it is working 
successfully we will probably need to go all the way to l-BFGS under 
double precision.

Let me point out that the molecules we are trying have longer than 1-4 
intramolecular interactions, and our [defaults ] section reads "1 2 yes 
0.5 0.8333", while the exclusion number in the [ moleculetype ] section 
is 3: as such, 1-4 interactions are scaled while longer ones should be 
fully taken into account.

As a first step, we made the energy optimization of a molecule, 
atoms+shells, without any restraint applied. So we just had to add the [ 
polarization ] section to the topology:

[ polarization ]
;  ai    sj  type    alpha (nm^3)
     1    10     1    500
     2    11     1    300
     3    12     1    500
(...)

We verified that the shells were moving related to "their" atoms 
because we made the calculations with different values of polarizability 
in the topology file (its magnitude was changed from e-3 to e+2 nm^3 for 
this testing purposes) and measured the shell distances respective to 
its atom. As expected, the distances between both runs were different 
(e-2 to e-3 nm) despite starting at the same initial positions, and as 
such we concluded that the minimization is acting over both atoms and 
shells due to only intramolecular interactions and in the absence of 
outer fields.

The second step was then to apply position restraints to all atoms, and 
*not* to any of the shells.

[ position_restraints ]
; ai  funct  fcx    fcy    fcz  restrains to a point
    1    1   10000000000  10000000000  10000000000
    2    1   10000000000  10000000000  10000000000
    3    1   10000000000  10000000000  10000000000
(...)
    9    1   10000000000  10000000000  10000000000

Then when we ran the energy optimization again, no matter the 
polarizability value applied (same range of magnitude, from e-3 to e+2 
nm^3), both atoms (as expected) and shells (reason of concern) did not 
move in absolute.

 From our point of view, applying the restraints to the atoms positions 
should not affect the shells that are (freely) bonded to them. However, 
we do not know if this was intentional from the developers for some 
reason, since we can understand that this is not a commonplace 
calculation.

We wish to know if this behavior is to be expected from the developers 
point of view (in which case, we will definitely move ahead towards 
studying and rewriting this piece of code by ourselves). If it is not, 
we will look forward to fill in a bug report (btw, 
http://bugzilla.gromacs.org/ seems to be offline right now, so we could 
not check if this is a known issue of some sort) and help each other.

Any ideas?

Thanks a lot in advance. Best Regards,

Jones