Minor issue regarding residue names.

[gduarter]

This is a relatively minor issue, but the name of the residue seems to affect the free energy calculation in unpredictable ways. In example1, all solvent molecules are labeled as 0C9 in the .yaml file and in SYSTEM.prmtop. The simulation runs normally but yields the strange result of 67.851 +- 0.128 kcal/mol.

example1.zip

In example2, all solvent molecules are labeled as XXX in the .yaml file and in SYSTEM.prmtop. A simulation run with these files yields the correct result of -5.140 +- 0.094 kcal/mol. The reference value for this solvation free energy is -5.2 +- 0.1 kcal/mol.

example2.zip

I think it would be interesting to either interrupt a simulation in the beginning if the residue name starts with a number or make Yank accept these unusual names. You can find more details in the attached files.

Thanks!

[jchodera]

@Lnaden : Can you look into this?

[jchodera]

Presumably, this is a LEaP issue, so examining the output of leap.log is likely necessary to sort out what's going on here. Examining the serialized system.xml (before alchemical modification) may also be helpful.

[jchodera]

@gduarter : If you still have them, could you ZIP up example setup/ directories from these and post them here?

[gduarter]

Here are the setup.log files for example1 and example2:
setup_1.zip
setup_2.zip

[jchodera]

Hm, I think we need the subdirectories too. We need the leap output logs for certain.

[andrrizzi]

I think I implemented a work-around in #819 that was released with YANK 0.19 to avoid Leap crashing when a number was in the residue name. If you're using a version >= 0.19 that might be the culprit.

[gduarter]

@jchodera: unfortunately there weren't any subdirectories in setup/, just setup.log. I don't know if it might help, but here's the everything but the trajectories:
example1_all.zip
example2_all.zip

@andrrizzi: I'm using 0.18.

[jchodera]

@andrrizzi : If you could get this fixed when you return, that would be great!

[andrrizzi]

This looks like an MDTraj problem. With

import mdtraj

topology = mdtraj.load_topology('SYSTEM.prmtop')
print('N atoms: ', topology.n_atoms)
print('N residues', topology.n_residues)
print('N ligand atoms: ', len(topology.select('resname LIG')))
print('N solvent atoms: ', len(topology.select('resname 0C9')))

I got

N atoms:  2822
N residues 157
N ligand atoms:  14
N solvent atoms:  0

I'll open an issue there.