michellab / BioSimSpace

Code and resources for the EPSRC BioSimSpace project.
https://biosimspace.org
GNU General Public License v3.0
77 stars 19 forks source link

engine="gromacs", gives strange gro file #196

Closed rongfengzou closed 3 years ago

rongfengzou commented 3 years ago

Dear BioSimSpace developer

I ran the workshops/advanced/alchemistry/answers/alchemistry_setup_solution.ipynb notebook. The outputs were saved to gromacs format. I viewed the gro file with VMD, but the gro file looks very strange. Moreover, gmx_mpi mdrun -deffnm gromacs gives errors. gromacs.gro.zip

Could you please help me deal with this issue?

Best regards,

lohedges commented 3 years ago

Hi there,

Could you provide a little more context please, i.e. which version of BioSimSpace are you using, what operating system, etc?

The GRO file will look strange when visualised since it contains the state of the merged system, i.e. it contains coordinates for atoms at both end points. If you want to visualise it in a sensible way, then you'd need to convert the merged molecule to its state at either lambda = 0 or lambda = 1, then visualise that.

Moreover, gmx_mpi mdrun -deffnm gromacs gives errors.

What were the errors? Could you post them here, or upload a zip file containing the working directory of the process?

rongfengzou commented 3 years ago

Dear BioSimSpace developer,

Here' re the code I used:

import BioSimSpace as BSS

lig1 = 'input1.mol2' lig2 = 'input2.mol2'

Load the protein and two ligands.

protein = BSS.IO.readMolecules("protein.pdb")[0] l1 = BSS.IO.readMolecules(lig1)[0] l2 = BSS.IO.readMolecules(lig2)[0]

gen parameters for protein

pro = BSS.Parameters.ff99SBildn(protein).getMolecule() ll1 = BSS.Parameters.gaff2(l1, net_charge=-1).getMolecule() ll2 = BSS.Parameters.gaff2(l2, net_charge=-1).getMolecule()

mapping = BSS.Align.matchAtoms(ll1, ll2)

ll1 = BSS.Align.rmsdAlign(ll1, ll2, mapping) merged = BSS.Align.merge(ll1, ll2, mapping, allow_ring_breaking=True,allow_ring_size_change=True ) system = merged + pro

Solvate in a 60 angstrom box of TIP3P water.

solvated = BSS.Solvent.tip3p(molecule=system, box=3[60BSS.Units.Length.angstrom], ion_conc=0.15)

Minimise the system.

minimised = BSS.Process.Gromacs(solvated, BSS.Protocol.Minimisation()).start().getSystem(block=True)

Equilibrate the system.

equilibrated = BSS.Process.Gromacs(minimised, BSS.Protocol.Equilibration()).start().getSystem(block=True)

Create the free energy protocol.

protocol = BSS.Protocol.FreeEnergy(runtime=5*BSS.Units.Time.nanosecond, num_lam=12)

Initialise the binding free energy object.

freenrg = BSS.FreeEnergy.Binding(solvated, protocol, work_dir="ti", engine='gromacs')

The version I used is 2019.3.0+58.ga1409d3.

Attached are the lambda_1.0000 and ligand mol2 files.

Best Regards, Rongfeng lambda_1.0000.tar.gz ligand.zip

rongfengzou commented 3 years ago

Dear BioSimSpace developer,

I think I solved the issue. It is OK to close it.

Best, Rongfeng

lohedges commented 3 years ago

Thanks for the update, I'll close this now.

lohedges commented 3 years ago

Just to note that you are using a very old version of BioSimSpace. If you're using conda, see the instructions in the README on the main page for details on how to ensure you get the most recent version.