Codecov Report

Merging #174 (35583cb) into develop (cc104f0) will decrease coverage by 0.01%.
The diff coverage is n/a.

@@             Coverage Diff             @@
##           develop     #174      +/-   ##
===========================================
- Coverage    78.90%   78.88%   -0.02%     
===========================================
  Files           11       11              
  Lines         1692     1691       -1     
  Branches       269      269              
===========================================
- Hits          1335     1334       -1     
- Misses         274      275       +1     
+ Partials        83       82       -1     

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

@orbeckst It only raised the error when gromacs wanted to replace some solvent molecules to add ions, which would recognize the molecule name from the atom section in tip files while the molecule name in gro or top files were based on the name defined in molecule section.

When dealing with neutral molecules (SAMPL), it wouldn't need to replace solvent molecules, so no error was raised. But when the molecule has charge, such as Rick's, the error is raised.

To test this we need a charged molecule. Maybe I add a sodium ion system for the tests?

Great that you added the tests. The octanol-ion ones fail for most GROMACS versions, though. The output is not very informative, though, so I have no idea what's going wrong.

I'm going to run a local test starting with gromacs wrapper's solvate function.

@orbeckst I think the problem is due to that gromacswrapper currently doesn't support np, and nn (number of negative and positive ions) input as well as using the concentration with a non-water solvent.

Here's the source code of gromacswrapper setup.py from line 451. We can see in the if concentration part, the selection inputs for generating an index file is only for water molecule.

    with in_dir(dirname):
        with open('none.mdp','w') as mdp:
            mdp.write('; empty mdp file\ninclude = {include!s}\nrcoulomb = 1\nrvdw = 1\nrlist = 1\n'.format(**mdp_kwargs))
        qtotgmx = cbook.grompp_qtot(f='none.mdp', o='topol.tpr', c=structure,
                                    p=topology, stdout=False, maxwarn=grompp_maxwarn)
        qtot = round(qtotgmx)
        logger.info("[{dirname!s}] After solvation: total charge qtot = {qtotgmx!r} = {qtot!r}".format(**vars()))

        if concentration != 0:
            logger.info("[{dirname!s}] Adding ions for c = {concentration:f} M...".format(**vars()))
            # target concentration of free ions c ==>
            #    N = N_water * c/c_water
            # add ions for concentration to the counter ions (counter ions are less free)
            #
            # get number of waters (count OW ... works for SPC*, TIP*P water models)
            rc,output,junk = gromacs.make_ndx(f='topol.tpr', o='ow.ndx',
                                              input=('keep 0', 'del 0', 'a OW*', 'name 0 OW', '', 'q'),
                                              stdout=False)
            groups = cbook.parse_ndxlist(output)
            gdict = {g['name']: g for g in groups}   # overkill...
            N_water = gdict['OW']['natoms']                  # ... but dict lookup is nice
            N_ions = int(N_water * concentration/CONC_WATER) # number of monovalents
        else:
            N_ions = 0

That makes sense. The concentration calculation wouldn't make any sense with non-water solvents. However, I am not even sure that the GW code works with all water models – do they all actually call the oxygen OW??

Does your change still work with the mixed solvent?

This issues looks as if it's tangled up with how to deal with ions properly. My understanding is that it does not limit the current functionality, right? Therefore, we can push this to a later milestone.

I think the current version of groamcs wrapper's solvate function with concentration input only works with some specific water models. For the test itself, we need a update for gromcas wrapper to support np, nn inputs, and we don't need to think about ion concentration in octanol.

It still works with met octanol. In the topology files of mixed octanol, water molecules were named as SOL which is the atomtype selected for being replaced with ions, and octanol was defined as OcOH``.

Can you please raise an issue for GW solvate to support nn/np and mention this PR, so that we don't forget? Thanks!

@VOD555 is this a PR we should get finished for a 0.8.0 release?

@orbeckst It's not necessary for normal use of MDPOW (for neutral molecules), so I think it can be finished after 0.8.0. And we still need a gromacswrapper supports nn and np input.