Failed to minimise perturbable system

[kexul]

Failure Information

'gmx grompp' reported the following errors:
ERROR 1 [file gromacs.top, line 393]:
  virtual site H (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.
  ERROR 2 [file gromacs.top, line 393]:
  virtual site H1 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.
  ERROR 3 [file gromacs.top, line 393]:
  virtual site H2 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.

Steps to Reproduce

import BioSimSpace as BSS

lig = BSS.IO.readPerturbableSystem('cor.rst7', 'top0.prm7', 'top1.prm7')

solvated_sys = BSS.Solvent.tip3p(lig, box=3*[8*BSS.Units.Length.nanometer], work_dir='solvation')

minimise_protocol = BSS.Protocol.Minimisation()
minimised_proc = BSS.Process.Gromacs(solvated_sys, minimise_protocol, work_dir='minimization')
minimised_proc.start()
minimised_proc.wait()
minimised_sys = minimised_proc.getSystem()

git2.zip

[lohedges]

Thanks. If possible, could you provide the original files and script used to create the merged ligand. I'd like to see if it works before writing the files to disk and reading back. That way I can check what's different between the two systems.

Cheers.

[lohedges]

Hmmm. As it says, it's because the dummy atoms have non-zero mass, i.e. it takes the mass from the end state with the heavier atom.

for atom in lig.getAtoms():
    print(atom._sire_object.property("element0"), atom._sire_object.property("mass0"))

Carbon (C, 6) 12.01 g mol-1
Nitrogen (N, 7) 14.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Nitrogen (N, 7) 14.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Nitrogen (N, 7) 14.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Nitrogen (N, 7) 14.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Bromine (Br, 35) 79.9 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Carbon (C, 6) 12.01 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
Hydrogen (H, 1) 1.008 g mol-1
dummy (Xx, 0) 1.008 g mol-1
dummy (Xx, 0) 1.008 g mol-1
dummy (Xx, 0) 1.008 g mol-1

I'm sure I've run simulations on the lambda=0 state using GROMACS before. When I get a chance, I'll check what's changed.

[lohedges]

Ah, I think I must have used a perturbation where there were no dummies at lambda=0. If the above is correct, then we wouldn't be able to simulate the lambda=0 system when dummies are present. Will see if I can think of a workaround. (It's possible that both end states could have dummies, so we can just choose the one without.)

[lohedges]

Okay, it's related to this issue. Essentially we need to flag the dummy as ptype A in order to simulate it with mass. I'll see if I can update this in BioSimSpace. If not, I'll need to adjust the Sire.IO.GroTop parser.

[lohedges]

Also, just to note, that having set the type to A, gmx grompp complains about the following:

NOTE 1 [file gromacs.top, line 393]:
  System has non-zero total charge: -0.000999
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.

It looks like the original ligands don't have integer charge. This should be done if they were parameterised by BioSimSpace. If not, you can call molecule._fixCharge() to do this for you.

[kexul]

Thanks for the quick response! I'm still finding my original input...

[lohedges]

Don't worry about that, it's not the source of the problem. (Should have mentioned this, whoops.)

[lohedges]

I've got it to work by not extracting the lambda=0 state prior to writing the GROMACS topology file. GROMACS just uses the lambda=0 state from the file when performing non-FEP simulations anyway, so there is no need. When writing the two-state system, Sire writes the correct ptype value for dummy atoms with mass, i.e. A.

The following should now work. (Note that I've fixed the charge of the ligands, but it's probably best to do this prior to merging.)

import BioSimSpace as BSS

lig = BSS.IO.readPerturbableSystem('cor.rst7', 'top0.prm7', 'top1.prm7')

# Extract the merged molecule and make sure the charge in each end state is
# integer valued.
merged = lig[0]
merged._fixCharge(property_map={"charge" : "charge0"})
merged._fixCharge(property_map={"charge" : "charge1"})
lig.updateMolecules(merged)

solvated_sys = BSS.Solvent.tip3p(lig, box=3*[8*BSS.Units.Length.nanometer], work_dir='solvation')

minimise_protocol = BSS.Protocol.Minimisation()
minimised_proc = BSS.Process.Gromacs(solvated_sys, minimise_protocol, work_dir='minimization')
minimised_proc.start()
minimised_proc.wait()
minimised_sys = minimised_proc.getSystem()

[kexul]

Thanks for the quick fix! It works very well.

I'm using Gromacs now as you suggest that it's more stable. But I still got nan error after minimized by Gromacs, if equilibrated by SOMD:

equilibration_protocol = BSS.Protocol.Equilibration()
equilibrated_proc = BSS.Process.Somd(minimised_sys, equilibration_protocol, work_dir='equilibration')
equilibrated_proc.start()
equilibrated_proc.wait()
equilibrated_sys = equilibrated_proc.getSystem()

In fact, I prefer SOMD since it's faster. Are there parameters I can tune to get rid of the nan issue?

[lohedges]

I'm not sure. It could be a case of needing to minimise for longer. (When I ran a test with your system the energy didn't converge to the desired tolerance after 10000 steps.) What happens if you re-minimise the minimised GROMACS system with SOMD before performing the equilibration? Does that blow up immediately too? It could be that SOMD/OpenMM can't handle the perturbed system for some reason, other than for a FEP simulation.

[lohedges]

Since SOMD/OpenMM load from AMBER format files, it could again be an issue with the unknown du type, which I assume is only used by SOMD internally for setting up the FEP simulations.

[kexul]

What happens if you re-minimise the minimised GROMACS system with SOMD before performing the equilibration? Does that blow up immediately too?

Thanks, I'll try it tomorrow.

It could be that SOMD/OpenMM can't handle the perturbed system for some reason, other than for a FEP simulation.

If that's the case,would it be supported in the future?😃

[lohedges]

If that's the case,would it be supported in the future?smiley

Well, there's no reason why it shouldn't be possible. I really hope that we can get it to work with all engines eventually.

For reference, SOMD can minimise the system (the output of the GROMACS minimisation), but still blows up when you try to equilibrate.

[kexul]

Well, there's no reason why it shouldn't be possible. I really hope that we can get it to work with all engines eventually.

Glad to hear that！

Thanks for the information😃

[kexul]

Hi @lohedges, I've done some tests in recent days, it seems like SOMD is somehow not compatible with the system minimized and equilibrated by Gromacs. NaN error was raised during the simulation of complex leg (free leg finished normally).

NaN or Inf has been generated along the simulation

Here is the code I've used to minimize and equilibrate the system, mostly copied from the tutorial:

import BioSimSpace as BSS
import argparse
import os
from pathlib import Path
import logging

logging.basicConfig(level=logging.INFO)

def get_box_length(sys, sys_type):
    box_min, box_max = sys.getAxisAlignedBoundingBox()
    box_size = [y - x for x, y in zip(box_min, box_max)]
    if sys_type == 'ligand':
        padding = 15 * BSS.Units.Length.angstrom
    elif sys_type == 'complex':
        padding = 8 * BSS.Units.Length.angstrom
    else:
        raise NotImplementedError(f'Not compatitable for this {sys_type}')
    box_length = max(box_size) + padding
    return box_length

def run_process(system, protocol, work_dir):
    """
    Run system with protocol, get the system afterward
    Parameters
    ----------
    system: sire system
    protocol: BioSimSpace Protocol

    Returns
    -------
    sire system
    """
    process = BSS.Process.Gromacs(system, protocol, work_dir=work_dir, ignore_warnings=True)
    process.addArgs({'-ntmpi':1, '-ntomp':'30', '-gpu_id':0})
    process.start()
    process.wait()
    system = process.getSystem()
    return system

def minimise_and_equilibrate(sys, leg='complex'):
    """
    Minimize the system and then equilibrate it
    Parameters
    ----------
    sys: sire system
    leg: complex or ligand

    Returns
    -------
    None
    """
    runtime_short_nvt = 50 # ps
    runtime_nvt = 200 # ps RESET to 50 after TESTING !
    runtime_npt = 500 # ps RESET to 200 after TESTING !

    logging.info(f'Solvating {leg}')
    solvated_sys = BSS.Solvent.tip3p(sys, box=3*[get_box_length(sys, leg)], work_dir=f'Solvation/{leg}')

    logging.info(f'Minimising {leg}')
    minimise_protocol = BSS.Protocol.Minimisation()
    minimised_sys = run_process(solvated_sys, minimise_protocol, work_dir=f'Minimization/{leg}')

    logging.info(f'Equilibrating {leg}')
    nvt_short = BSS.Protocol.Equilibration(runtime=runtime_short_nvt*BSS.Units.Time.picosecond,
                                           temperature_start=0*BSS.Units.Temperature.kelvin,
                                           temperature_end=300*BSS.Units.Temperature.kelvin,
                                           restraint="all")
    equil1 = run_process(minimised_sys, nvt_short, work_dir=f'Nvt_short/{leg}')

    if leg == 'complex':
        nvt_backbone = BSS.Protocol.Equilibration(
            runtime=runtime_nvt*BSS.Units.Time.picosecond,
            temperature=300*BSS.Units.Temperature.kelvin,
            restraint="backbone"
        )
        equil1 = run_process(equil1, nvt_backbone, work_dir=f'Nvt_backbone/{leg}')

    nvt = BSS.Protocol.Equilibration(
        runtime=runtime_nvt*BSS.Units.Time.picosecond,
        temperature=300*BSS.Units.Temperature.kelvin,
    )
    equil2 = run_process(equil1, nvt, work_dir=f'Nvt/{leg}')

    npt_heavy = BSS.Protocol.Equilibration(
        runtime=runtime_npt*BSS.Units.Time.picosecond,
        pressure=1*BSS.Units.Pressure.atm,
        temperature=300*BSS.Units.Temperature.kelvin,
        restraint="heavy",
    )
    equil3 = run_process(equil2, npt_heavy, work_dir=f'Npt_heavy/{leg}')

    npt = BSS.Protocol.Equilibration(
        runtime=runtime_npt*BSS.Units.Time.picosecond,
        pressure=1*BSS.Units.Pressure.atm,
        temperature=300*BSS.Units.Temperature.kelvin,
    )
    equilibrated_sys = run_process(equil3, npt, work_dir=f'Npt/{leg}')

    logging.info(f'Creating SOMD binding protocal for perturbation')
    protocol = BSS.Protocol.FreeEnergy(num_lam=32, report_interval=100000, restart_interval=50000)
    if leg == 'complex':
        BSS.FreeEnergy.Relative(equilibrated_sys, protocol, engine='SOMD',
                                work_dir='bound', setup_only=True)
    elif leg == 'ligand':
        BSS.FreeEnergy.Relative(equilibrated_sys, protocol, engine='SOMD',
                                work_dir='free', setup_only=True)

def create_system(working_dir):
    os.chdir(working_dir)
    complex_template = BSS.IO.readMolecules('complex_template.*')
    rst = str(list(Path('.').glob('*.rst7'))[0])
    p0 = str(list(Path('.').glob('*0.prm7'))[0])
    p1 = str(list(Path('.').glob('*1.prm7'))[0])
    morph = BSS.IO.readPerturbableSystem(rst, p0, p1)

    minimise_and_equilibrate(morph, leg='ligand')
    minimise_and_equilibrate(complex_template+morph, leg='complex')

if __name__ == '__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument("-w", "--workdir", type=str, help="Path to where "
                                                          "simulation files are located")
    args = parser.parse_args()
    create_system(args.workdir)
# then we run somd-freenrg in the workdir/bound and workdir/free. 

I've confirmed that the simulation would end without any error if the system were prepared by directly fed to BSS.FreeEnergy.Relative, that is :

import BioSimSpace as BSS
import argparse
import os
from pathlib import Path
import logging

logging.basicConfig(level=logging.INFO)

def create_system(working_dir):
    os.chdir(working_dir)
    complex_template = BSS.IO.readMolecules('complex_template.*')
    rst = str(list(Path('.').glob('*.rst7'))[0])
    p0 = str(list(Path('.').glob('*0.prm7'))[0])
    p1 = str(list(Path('.').glob('*1.prm7'))[0])
    morph = BSS.IO.readPerturbableSystem(rst, p0, p1)

    logging.info(f'Creating SOMD binding protocal for perturbation')
    protocol = BSS.Protocol.FreeEnergy(num_lam=32, report_interval=100000, restart_interval=50000)

    solvated_complex = BSS.Solvent.tip3p(complex_template+morph, box=3*[get_box_length(complex_template+morph, 'complex')], work_dir=f'Solvation/complex')
    solvated_ligand = BSS.Solvent.tip3p(morph, box=3*[get_box_length(morph, 'ligand')], work_dir=f'Solvation/ligand')

    BSS.FreeEnergy.Relative(solvated_complex, protocol, engine='SOMD',
                            work_dir='bound', setup_only=True)
    BSS.FreeEnergy.Relative(solvated_ligand, protocol, engine='SOMD',
                        work_dir='free', setup_only=True)

if __name__ == '__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument("-w", "--workdir", type=str, help="Path to where "
                                                          "simulation files are located")
    args = parser.parse_args()
    create_system(args.workdir)
# then we run somd-freenrg in the workdir/bound and workdir/free. 

Here are the input files: git.zip

[kexul]

Two new findings:

1. If I re-run the simulation manually after it threw NaN, it may end normally.
2. In the log file of SOMD, the energy before minimize is quite large, maybe Gromacs minimization is not working.

   ###=======================Minimisation========================###
   Running minimisation.
   Energy before the minimisation: 2.52998e+09 kcal mol-1
   Tolerance for minimisation: 1
   Maximum number of minimisation iterations: 1000
   Energy after the minimization: -125968 kcal mol-1
   Energy minimization done.
   ###===========================================================###

[lohedges]

Thanks @kexul. I'm wondering if GROMACS actually isn't simulating the lambda=0 system, rather some intermediate state. Could you extract the "equlibrated" complex from the GROMACS simulations and visualise it to see if anything looks funky, e.g. bad contacts.

I'll also try running the GROMACS minimisation/equilibration in the old way, i.e. extracting just the lambda=0 state, and fixing the ptype in the topology file.

One other thought is the way in which SOMD/GROMACS handle potential terms for the perturbable system. Perhaps there is something that is "soft" in the GROMACS potential that is being applied in a "hard" way with SOMD. @jmichel80: Do you have any ideas here?

[kexul]

Could you extract the "equlibrated" complex from the GROMACS simulations and visualise it to see if anything looks funky, e.g. bad contacts.

Wow! Thanks for your suggestion @lohedges , the equilibrated complex is quite weird.

[lohedges]

Some comments related to GROMACS / SOMD differences can be found in this thread. It sounds like things should be the same at lambda=0.

[jmichel80]

Hi @kexul we occasionally see NaN during run time for large perturbations (e.g ring growht) although code changes in SOMD done last year helped reduce a lot occurences. The single point energy before minimisation can be high for lambda > 0 if atoms starting as dummies are clashing with solvent or protein atoms (the initial positioning of dummy atoms by BSS doesn't care about the rest of the system). If the system was equilibrated with a Hamiltonian equivalent to the lambda 0 state (apart from additional dummy atom bonded interactions) the single point energy should already be negative. It will go down a bit because the energy minimisation will quenched a thermally equilibrated system, but it shouldn't start with very positive values (which normally indicate repulsive clashes from a LJ component).

For the visulalisation, are you using an orthorhombic box or a triclinic one ? Gromacs snapshots have to be imaged carefully to avoid PBC visualisation artefacts.

[lohedges]

Okay, I modified the Sire.IO.GroTop parser so that it writes ptype A for dummy atoms if a molecule was formerly peturbable, i.e. we can convert a perturbable molecule to a lambda end state using molecule._toRegularMolecule() and write that single state to file and run it with GROMACS.

I then performed single-point energy calculations to compare the energy of the topology containing one state to the one containing both states, which is what is written at present. The energies agree perfectly so it appears that, for a non-FEP simulation, GROMACS is indeed simulating the lambba=0 state when you give it a topology file containing both states. (Which is what I thought was the case)

[kexul]

For the visulalisation, are you using an orthorhombic box or a triclinic one ? Gromacs snapshots have to be imaged carefully to avoid PBC visualisation artefacts.

Hi @jmichel80, I used the code below to solvate the system, I guess the box should be orthorhombic.

solvated_sys = BSS.Solvent.tip3p(sys, box=3*[get_box_length(sys, leg)], work_dir=f'Solvation/{leg}')

@lohedges, I found that the system became a mess after Minimization.

Here are all of the intermediate files during the Minimization and Equilibration procedure. https://drive.google.com/file/d/1f71EAHquiGikOHalAdNx_cJVCKwIRLr_/view?usp=sharing

[kexul]

I tried printing the potential energy during Gromacs Process, that is:

def run_process(system, protocol, work_dir):
    """
    Run system with protocol, get the system afterwards
    Parameters
    ----------
    system: sire system
    protocol: BioSimSpace Protocol

    Returns
    -------
    sire system
    """
    process = BSS.Process.Gromacs(system, protocol, work_dir=work_dir, ignore_warnings=True)
    logging.info(f'Before: {process.getCurrentPotentialEnergy()}')
    process.addArgs({'-ntmpi':1, '-ntomp':'30', '-gpu_id':0})
    process.start()
    process.wait()
    logging.info(f'After: {process.getCurrentPotentialEnergy()}')
    system = process.getSystem()
    return system

It showed:

INFO:root:Solvating ligand

INFO:root:Minimising ligand

INFO:root:Before: None

INFO:root:After: -8922.8728 kcal/mol

INFO:root:Equilibrating ligand

INFO:root:Before: None

INFO:root:After: -7415.4876 kcal/mol

INFO:root:Before: None

INFO:root:After: -7221.1281 kcal/mol

INFO:root:Before: None

INFO:root:After: -7203.0832 kcal/mol

INFO:root:Before: None

INFO:root:After: -7137.3088 kcal/mol

INFO:root:Creating SOMD binding protocal for perturbation

Gromacs seems to be minimizing the system correctly. But when I run somd-freenrg in the lambda_0.0000 folder, it showed that :

###================Setting up calculation=====================###
New run. Loading input and creating restart
lambda is 0.0
Create the System...
Selecting dummy groups
Creating force fields...
Setting up the simulation with random seed 580536
Setting up moves...
Created a MD move that uses OpenMM for all molecules on 0
Generated random seed number 580536
Saving restart
Setting up sim file.
There are 2244 atoms in the group
###===========================================================###

###=======================Minimisation========================###
Running minimisation.
Energy before the minimisation: 4.8221e+07 kcal mol-1
Tolerance for minimisation: 1
Maximum number of minimisation iterations: 1000

==============================================================
Sending anonymous Sire usage statistics to http://siremol.org.
For more information, see http://siremol.org/analytics
To disable, set the environment variable 'SIRE_DONT_PHONEHOME' to 1
To see the information sent, set the environment variable
SIRE_VERBOSE_PHONEHOME equal to 1. To silence this message, set
the environment variable SIRE_SILENT_PHONEHOME to 1.
==============================================================

Energy after the minimization: -8992.34 kcal mol-1
Energy minimization done.
###===========================================================###

Something might go wrong during the system conversion.

[lohedges]

Thanks, I'll run some single-point calculations here to see if I can figure out what's going on. Just a note that the potential energy will always be None before running a simulation, since this is extracted from the simulation log file, rather than computed from the system itself.

[lohedges]

I've taken your morph and run single-point energy calculation in vacuum using AMBER, GROMACS, and OpenMM. (For AMBER I converted the du atom type to EP.) The energies agree within tolerance for AMBER and GROMACS (bonded terms and potential energy) but the potential energy disagrees for OpenMM. I'm not sure if it calculates things in a different way, so I'm checking for a regular, non perturbable, molecule in vacuum too. Note that AMBER and OpenMM read the exact same AMBER input files, so there shouldn't be an issue with the GROMACS/AMBER conversion if both GROMACS and AMBER are giving the same single-point energies. Note also that I've excluded solvated calculations due to the different way in which AMBER and GROMACS treat fast three-point water molecules, i.e. explicit bonds vs restraint.

[kexul]

Hi, @lohedges, would you mind share your script used to calculate single-point energy? I would like to implement a PertrubableSystem reader for Gromacs.

[lohedges]

Sure thing. I normally do something like:

import BIoSimSpace as BSS

# Set up a vacuum system of some sort.
system = ...

# Create a minimisation protocol with a single step.
protocol = BSS.Protocol.Minimisation(steps=1)

# Create a process to run this with AMBER.
amb_process = BSS.Process.Amber(system, protocol, work_dir="amb_test")

# Create a process to run this with GROMACS.
gmx_process = BSS.Process.Gromacs(system, protocol, work_dir="gmx_test")

# To get a better comparison, set the number of steps to zero for GROMACS. We don't support this
# in the Protocol by default, since most engines don't allow zero steps.
config = gmx_process.getConfig()
config[1] = "nsteps = 0"
gmx_process.setConfig(config)

# Now run the minimisations.
amb_process.start()
amb_process.wait()
gmx_process.start()
gmx_process.wait()

# Now print out the energy terms that you want to check. Note that different engines
# output different information for different protocols, so you might not be able to compare
# directly through BioSimSpace. Also, GROMACS splits dihedral enegies into proper and
# impropers, whereas AMBER doesn't. I would suggest reading the values from
# amb_test/amber.nrg and gmx_test/gromacs.log and comparing manually.
print(amb_process.getBondEnergy(), gmx_process.getBondEnergy().kcal_per_mol())

Note that reading GROMACS format perturbable topologies wouldn't allow us to run a simulation directly with SOMD, since it requires some information that is specific to the AMBER format. That said, most of this could probably be recovered with a bit of work.

[lohedges]

Our GROMACS topology parser is here. Currently this can only write perturbable topologies, whereas it can read/write all other topologies. I had always planned to add a reader too, but it is much more complicated than writing.

[kexul]

Note that reading GROMACS format perturbable topologies wouldn't allow us to run a simulation directly with SOMD, since it requires some information that is specific to the AMBER format. That said, most of this could probably be recovered with a bit of work.

Thanks for pointing that out！😃

[kexul]

Hi @lohedges , I found that I'm getting this error again with the same input with the newest code.

'gmx grompp' reported the following errors:
ERROR 1 [file gromacs.top, line 393]:
  virtual site H (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.
  ERROR 2 [file gromacs.top, line 393]:
  virtual site H1 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.
  ERROR 3 [file gromacs.top, line 393]:
  virtual site H2 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.

[lohedges]

Could you check the pytpe for the atoms in your topology file. They should all read A. (You could also updload it here if needed.) Looking at the Anaconda page it doesn't look like any of the recent packages have been downloaded, so perhaps you are still using an outdated version. (Ignore this if you are using a source install, or the binary.)

[lohedges]

It might be that you need to update Sire, not BioSimSpace, since it is the Sire GroTop parser that was updated. (Apologies for not making this obvious earlier.)

[kexul]

Ignore this if you are using a source install, or the binary

Yes, I just cloned the newest code from github, not from Anaconda. Thanks, I'll try to update sire then.

[kexul]

virtual site H2 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008

After updating sire, this error no longer exists.

One thing to note, the latest version of biosimspace seems to conflict with the latest version of sire, here is the output of installing them:

conda create -n biosimspace -c conda-forge -c omnia -c michellab biosimspace=2020.1.0 sire=2021.1.0                                                                                          

Collecting package metadata (current_repodata.json): done                                                       
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.     
Collecting package metadata (repodata.json): done                                                               
Solving environment: -                                                                                          
Found conflicts! Looking for incompatible packages.                                                             
This can take several minutes.  Press CTRL-C to abort.                                                          
failed                                                                                                          

UnsatisfiableError: The following specifications were found to be incompatible with each other:                 

Output in format: Requested package -> Available versions                                                       

Package sire conflicts for:                                                                                     
biosimspace=2020.1.0 -> sire[version='>=2020.1.0,<2021.0a0']                                                    
sire=2021.1.0The following specifications were found to be incompatible with your system:                       

  - feature:/linux-64::__glibc==2.17=0                                                                          
  - feature:|@/linux-64::__glibc==2.17=0                                                                        
  - biosimspace=2020.1.0 -> libgcc-ng[version='>=7.5.0'] -> __glibc[version='>=2.17']                           

Your installed version is: 2.17       

[kexul]

Getting this error again when I use readMolecules() to construct the system, here is the script:

import BioSimSpace as BSS

lig = BSS.IO.readPerturbableSystem('cor.rst7', 'top0.prm7', 'top1.prm7') 
BSS.Solvent.tip3p(lig, box=3*[8*BSS.Units.Length.nanometer], work_dir='solvation')
solvated_sys = BSS.IO.readMolecules(['solvation/solvated_ions.gro', 'solvation/solvated.top'])
minimise_protocol = BSS.Protocol.Minimisation()
minimised_proc = BSS.Process.Gromacs(solvated_sys, minimise_protocol, work_dir='minimization')
minimised_proc.start()
minimised_proc.wait()
minimised_sys = minimised_proc.getSystem()

Here is the output:

'gmx grompp' reported the following errors:
ERROR 1 [file gromacs.top, line 378]:
  virtual site H (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.
  ERROR 2 [file gromacs.top, line 378]:
  virtual site H1 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.
  ERROR 3 [file gromacs.top, line 378]:
  virtual site H2 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008
  (state B)
  Check your topology.

If the return value of solvation was passed to process directly, everything is fine:

import BioSimSpace as BSS

lig = BSS.IO.readPerturbableSystem('cor.rst7', 'top0.prm7', 'top1.prm7') 
solvated_sys = BSS.Solvent.tip3p(lig, box=3*[8*BSS.Units.Length.nanometer], work_dir='solvation')
minimise_protocol = BSS.Protocol.Minimisation()
minimised_proc = BSS.Process.Gromacs(solvated_sys, minimise_protocol, work_dir='minimization', ignore_warnings=True)
minimised_proc.start()
minimised_proc.wait()
minimised_sys = minimised_proc.getSystem()

[lohedges]

Why do you need to do the first version? If you use BioSimSpace directly, as you are in the second case, the waters are copied back into the original system, which contains all of the properties required to identify a perturbable system. Reading the files from the working directory loses all of this information, so the GroTop parser doesn't know to treat perturbable dummy atoms differently, i.e. to use ptype=A.

[kexul]

Why do you need to do the first version? If you use BioSimSpace directly, as you are in the second case, the waters are copied back into the original system, which contains all of the properties required to identify a perturbable system. Reading the files from the working directory loses all of this information, so the GroTop parser doesn't know to treat perturbable dummy atoms differently, i.e. to use ptype=A.

Oops！I forgot that's a perturbable sysytem 😂

[kexul]

virtual site H2 (Res LIG-1) has non-zero mass 1.008 (state A) / 1.008

After updating sire, this error no longer exists.

One thing to note, the latest version of biosimspace seems to conflict with the latest version of sire, here is the output of installing them:

conda create -n biosimspace -c conda-forge -c omnia -c michellab biosimspace=2020.1.0 sire=2021.1.0                                                                                          

Collecting package metadata (current_repodata.json): done                                                       
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.     
Collecting package metadata (repodata.json): done                                                               
Solving environment: -                                                                                          
Found conflicts! Looking for incompatible packages.                                                             
This can take several minutes.  Press CTRL-C to abort.                                                          
failed                                                                                                          

UnsatisfiableError: The following specifications were found to be incompatible with each other:                 

Output in format: Requested package -> Available versions                                                       

Package sire conflicts for:                                                                                     
biosimspace=2020.1.0 -> sire[version='>=2020.1.0,<2021.0a0']                                                    
sire=2021.1.0The following specifications were found to be incompatible with your system:                       

  - feature:/linux-64::__glibc==2.17=0                                                                          
  - feature:|@/linux-64::__glibc==2.17=0                                                                        
  - biosimspace=2020.1.0 -> libgcc-ng[version='>=7.5.0'] -> __glibc[version='>=2.17']                           

Your installed version is: 2.17       

Hi @lohedges , would you mind fixing this?

[lohedges]

I'm not sure how that's possible, since they were both built using the exact same conda environment. Will check when I'm in front of the computer.

[lohedges]

It's also complaining about something that is apparently satisfied by the information in it's error report.

[lohedges]

It works fine for me, both with conda and mamba. I'm using Miniconda, which has conda version 4.9.2:

/opt/miniconda3/bin/conda --version
conda 4.9.2

[lohedges]

For reference, I used the following command, i.e. not pinning versions:

conda create -n test -c conda-forge -c omnia -c michellab/label/dev biosimspace 

(I checked the package plan contained the most recent versions of Sire and BioSimSpace, though.)

[lohedges]

Ah, looking at your command you are not using the development label for the michellab channel. The BioSImSpace package on the main label is older, hence the reason for the incompatibility. (The latest devel version of BioSimSpace has always been built against the latest devel version of Sire.) However, it should still work (as long as the Sire API is compatible, which it is) and it also doesn't explain the completely unhelpful (and incorrect) error message.

[kexul]

Thank you @lohedges, dev label worked for me too!🍻