Fix oct topology molecule name

[VOD555]

Fix #154
In itp files for dry octanol, the molecule name in atom section doesn't match the one in molecule section.
So, here, we replace OcOH with SOL to match the molecule names.

[codecov]

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     

Impacted Files	Coverage Δ
mdpow/equil.py	80.07% <0.00%> (-0.39%)	⬇️
mdpow/fep.py	69.18% <0.00%> (+0.12%)	⬆️

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

[orbeckst]

Why did this work in the past? Did we use custom top files??

Can you add tests, please? We probably only need to run grompp, something like the first stage of the equilibrium workflow.

[VOD555]

@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?

[orbeckst]

There are various issues related to ions #97 #85 and in the docs I said that charged molecules are currently no supported. Let's not start with ions right here, at least not as a test case for a FEP calculation.

For a test it then you could add a test that builds an octanol system (using Goct.solvate() etc) with options to add 4 Na and 4 Cl — we don't have to run the system, just checking that it generates it and that it contains the right number of ions afterwards.

Please also update CHANGES.

[orbeckst]

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.

[VOD555]

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

[VOD555]

@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

[orbeckst]

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?

MDPOW/mdpow/equil.py

Lines 637 to 649 in 43485dd

	def _setup_solvate(self, **kwargs):
	sol = gromacs.setup.solvate_sol(**kwargs)
	with in_dir(self.dirs.solvation, create=False):
	u = mda.Universe('solvated.gro')
	octanol = u.select_atoms('resname OcOH')
	n = octanol.n_residues
	with in_dir(self.dirs.topology, create=False):
	gromacs.cbook.edit_txt(self.files.topology,
	[('OcOH 1', '1', n)])
	ionkwargs = kwargs
	ionkwargs['struct'] = sol['struct']
	params = gromacs.setup.solvate_ion(**ionkwargs)
	return params

seems to say that the mixed-solvent code expects OcOH.

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.

[VOD555]

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``.

[orbeckst]

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

[orbeckst]

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

[VOD555]

@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.