search

choderalab / yank

An open, extensible Python framework for GPU-accelerated alchemical free energy calculations.
http://getyank.org
MIT License
176 stars 70 forks source link

Error with !Combinatorial usage with epik protonation states #549

Closed jchodera closed 7 years ago

jchodera commented 7 years ago

I tried to specify a !Combinatorial tag to set up several protonation states in the epik: block in a variant of the old abl-imatinib-explicit example:


---
options:
  minimize: yes
  verbose: yes
  output_dir: .
  number_of_iterations: 100
  nsteps_per_iteration: 500
  temperature: 300*kelvin
  pressure: 1*atmosphere
  restraint_type: Boresch
  softcore_beta: 0.0
  timestep: 2.0*femtoseconds
  platform: OpenCL

molecules:
  Abl:
    filepath: input/2HYY-pdbfixer.pdb
  STI:
    filepath: input/STI02.mol2
    epik:
      select: !Combinatorial [0, 1, 2]
      ph: 7.4
      ph_tolerance: 0.7
      tautomerize: no
    openeye:
      quacpac: am1-bcc
    antechamber:
      charge_method: null

solvents:
  pme:
    nonbonded_method: PME
    nonbonded_cutoff: 9*angstroms
    clearance: 9*angstroms
    positive_ion: Na+
    negative_ion: Cl-

systems:
  Abl-STI:
    receptor: Abl
    ligand: STI
    solvent: pme
    leap:
      parameters: [leaprc.ff14SB, leaprc.gaff, frcmod.ionsjc_tip3p]

protocols:
  absolute-binding:
    complex:
      alchemical_path:
        lambda_restraints:     [0.00, 0.25, 0.50, 0.75, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.000, 1.00, 1.000, 1.00, 1.00]
        lambda_electrostatics: [1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.000, 0.00, 0.000, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.075, 0.05, 0.025, 0.01, 0.00]
    solvent:
      alchemical_path:
        lambda_electrostatics: [1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.00]

experiments:
  system: Abl-STI
  protocol: absolute-binding

The resulting error was:

[chodera@mskcc-ln1 abl-imatinib]$ cat ~/abl-imatinib-explicit.o7782194
Warning: no access to tty (Bad file descriptor).
Thus no job control in this shell.
Running simulation via MPI...
2016-11-09 01:25:20,615: MPI enabled.
2016-11-09 01:25:20,642: Setting up the systems for Abl, STI_2 and pme
Invalid Range specified. Use '-h' option for usage
FATAL m2io_file_new(): Can't get first token from file: /cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.mae. File is empty.
m2io_open failed for file /cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.mae.
2016-11-09 01:25:37,341: ERROR - openmoltools.schrodinger - b'ERROR: ERROR: output from the converter is missing: /cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.sdf\n'
2016-11-09 01:25:37,341: ERROR - openmoltools.schrodinger - Command '['/cbio/jclab/share/schrodinger/schrodinger2016-3/utilities/structconvert', '-imae', '/cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.mae', '-osd', '/cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.sdf']' returned non-zero exit status 1

Am I not allowed to specify a !Combinatorial option there, or does this have to appear in a specific order?

jchodera commented 7 years ago

I tried a single protonation state per molecule:

molecules:                                                                                                                                                                                                                                                                                                                                                                
  Abl:                                                                                                                                                                                                                                                                                                                                                                    
    filepath: input/2HYY-pdbfixer.pdb                                                                                                                                                                                                                                                                                                                                     
  STI0:                                                                                                                                                                                                                                                                                                                                                                   
    filepath: input/STI02.mol2                                                                                                                                                                                                                                                                                                                                            
    epik:                                                                                                                                                                                                                                                                                                                                                                 
      select: 0                                                                                                                                                                                                                                                                                                                                                           
      ph: 7.4                                                                                                                                                                                                                                                                                                                                                             
      ph_tolerance: 0.7                                                                                                                                                                                                                                                                                                                                                   
      tautomerize: no                                                                                                                                                                                                                                                                                                                                                     
    openeye:                                                                                                                                                                                                                                                                                                                                                              
      quacpac: am1-bcc                                                                                                                                                                                                                                                                                                                                                    
    antechamber:                                                                                                                                                                                                                                                                                                                                                          
      charge_method: null

but it looks like either the molecule names cannot contain numbers or the molecule is loaded into a unit named after the mol2 substructure name. From the failure solvent.leap.log:

> # Load molecules
> STI0 = loadMol2 moli1.mol2
Loading Mol2 file: ./moli1.mol2
Reading MOLECULE named MOL
> 
> # Solvate systems
> addIons2 STI0 Na+ 0
addIons: Argument #2 is type String must be of type: [unit]
jchodera commented 7 years ago

Wait, something else must be wrong.

jchodera commented 7 years ago

Nope, I can't even get the basic YAML to work without epik:

---
options:
  minimize: yes
  verbose: yes
  output_dir: .
  number_of_iterations: 100
  nsteps_per_iteration: 500
  temperature: 300*kelvin
  pressure: 1*atmosphere
  restraint_type: Boresch
  softcore_beta: 0.0
  timestep: 2.0*femtoseconds
  platform: OpenCL

molecules:
  Abl:
    filepath: input/2HYY-pdbfixer.pdb
  STI:
    filepath: input/STI02.mol2
    openeye:
      quacpac: am1-bcc
    antechamber:
      charge_method: null

solvents:
  pme:
    nonbonded_method: PME
    nonbonded_cutoff: 9*angstroms
    clearance: 9*angstroms
    positive_ion: Na+
    negative_ion: Cl-

systems:
  Abl-STI:
    receptor: Abl
    ligand: STI
    solvent: pme
    leap:
      parameters: [leaprc.ff14SB, leaprc.gaff, frcmod.ionsjc_tip3p]

protocols:
  absolute-binding:
    complex:
      alchemical_path:
        lambda_restraints:     [0.00, 0.25, 0.50, 0.75, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.000, 1.00, 1.000, 1.00, 1.00]
        lambda_electrostatics: [1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.000, 0.00, 0.000, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.075, 0.05, 0.025, 0.01, 0.00]
    solvent:
      alchemical_path:
        lambda_electrostatics: [1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.00]

experiments:
  system: Abl-STI
  protocol: absolute-binding
Lnaden commented 7 years ago

leaprc.ff14SB got moved to oldff/leaprc.ff14SB, It probably cant find the ion names because they were not loaded correctly

jchodera commented 7 years ago

Thanks! That fixed one problem when I don't include the epik block, but when I add it:

---
options:
  minimize: yes
  verbose: yes
  output_dir: .
  number_of_iterations: 100
  nsteps_per_iteration: 500
  temperature: 300*kelvin
  pressure: 1*atmosphere
  restraint_type: Boresch
  softcore_beta: 0.0
  timestep: 2.0*femtoseconds
  platform: OpenCL

molecules:
  Abl:
    filepath: input/2HYY-pdbfixer.pdb
  STI:
    filepath: input/STI02.mol2
    epik:
      select: !Combinatorial [1, 2, 3]
      ph: 7.4
      ph_tolerance: 0.7
      tautomerize: no
    openeye:
      quacpac: am1-bcc
    antechamber:
      charge_method: null

solvents:
  pme:
    nonbonded_method: PME
    nonbonded_cutoff: 9*angstroms
    clearance: 9*angstroms
    positive_ion: Na+
    negative_ion: Cl-

systems:
  Abl-STI:
    receptor: Abl
    ligand: STI
    solvent: pme
    leap:
      parameters: [oldff/leaprc.ff14SB, leaprc.gaff, frcmod.ionsjc_tip3p]

protocols:
  absolute-binding:
    complex:
      alchemical_path:
        lambda_restraints:     [0.00, 0.25, 0.50, 0.75, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.000, 1.00, 1.000, 1.00, 1.00]
        lambda_electrostatics: [1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.000, 0.00, 0.000, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.075, 0.05, 0.025, 0.01, 0.00]
    solvent:
      alchemical_path:
        lambda_electrostatics: [1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.00]

experiments:
  system: Abl-STI
  protocol: absolute-binding

I end up with 

[chodera@gpu-1-17 abl-imatinib]$ yank script --yaml=yank.yaml
2016-11-09 09:04:46,420: MPI disabled.
2016-11-09 09:04:46,445: Setting up the systems for Abl, STI_1 and pme
2016-11-09 09:05:19,894: 4-[(4-methylpiperazin-4-ium-1-yl)methyl]-~{N}-[4-methyl-3-[[4-(3-pyridyl)pyrimidin-2-yl]amino]phenyl]benzamide
2016-11-09 09:06:11,474: Setting up the systems for Abl, STI_2 and pme
Invalid Range specified. Use '-h' option for usage
FATAL m2io_file_new(): Can't get first token from file: /cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.mae. File is empty.
m2io_open failed for file /cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.mae.
2016-11-09 09:06:18,328: ERROR - openmoltools.schrodinger - b'ERROR: ERROR: output from the converter is missing: /cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.sdf\n'
2016-11-09 09:06:18,328: ERROR - openmoltools.schrodinger - Command '['/cbio/jclab/share/schrodinger/schrodinger2016-3/utilities/structconvert', '-imae', '/cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.mae', '-osd', '/cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.sdf']' returned non-zero exit status 1
Traceback (most recent call last):
  File "/cbio/jclab/home/chodera/miniconda3/bin/yank", line 9, in <module>
    load_entry_point('yank==0.12.0', 'console_scripts', 'yank')()
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/yank/cli.py", line 74, in main
    dispatched = getattr(commands, command).dispatch(command_args)
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/yank/commands/script.py", line 59, in dispatch
    yaml_builder.build_experiments()
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/yank/yamlbuild.py", line 1282, in build_experiments
    self._setup_experiments()
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/yank/yamlbuild.py", line 2011, in _setup_experiments
    self._db.get_system(system_id)
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/yank/yamlbuild.py", line 656, in get_system
    0, system_parameters, solvent_id, ligand_id)
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/yank/yamlbuild.py", line 964, in _setup_system
    self._setup_molecules(*molecule_ids)  # Be sure molecules are set up
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/yank/yamlbuild.py", line 866, in _setup_molecules
    omt.schrodinger.run_structconvert(epik_mae_file, epik_sdf_file)
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/openmoltools/schrodinger.py", line 74, in _need_schrodinger
    return func(*args, **kwargs)
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/openmoltools/schrodinger.py", line 164, in run_structconvert
    run_and_log_error(cmd)
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/openmoltools/schrodinger.py", line 43, in run_and_log_error
    raise e
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/site-packages/openmoltools/schrodinger.py", line 39, in run_and_log_error
    output = subprocess.check_output(command)
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/subprocess.py", line 626, in check_output
    **kwargs).stdout
  File "/cbio/jclab/home/chodera/miniconda3/lib/python3.5/subprocess.py", line 708, in run
    output=stdout, stderr=stderr)
subprocess.CalledProcessError: Command '['/cbio/jclab/share/schrodinger/schrodinger2016-3/utilities/structconvert', '-imae', '/cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.mae', '-osd', '/cbio/jclab/home/chodera/github/choderalab/yank/yank-examples/examples/binding/abl-imatinib/setup/molecules/STI_2/STI_2-epik.sdf']' returned non-zero exit status 1
andrrizzi commented 7 years ago

Can you take a look at the epik log? The most likely explanation to me is that epik does not generate 4 different protonation states (select is 0-based). This

Invalid Range specified. Use '-h' option for usage

could be the error message of maesubset, which is used to extract the particular structure you index.

jchodera commented 7 years ago

I'm not quite sure how to find the epik log. Maybe we should see if epik returns an error code on exit and then dump the epik log to the terminal if that happens?

andrrizzi commented 7 years ago

There should be an epik.log file (or some similar name) in the setup/molecules/STIx/ folder.

jchodera commented 7 years ago

Aha! Thanks!

jchodera commented 7 years ago

OK, the issue is that epik only gives two protonation states for a pH tolerance of 0.7.

Interestingly, epik is run multiple times---once for each set combinatorial option. Since it's so fast, this doesn't impact performance. But we should probably increase the pH tolerance in our example documentation by a lot so people don't run into this problem. It shouldn't slow things down any.

andrrizzi commented 7 years ago

Interestingly, epik is run multiple times

Yep! That's expected with !Combinatorial experiments. Effectively that keyword generate 3 different molecules and the pipeline doesn't know that the Epik result will be the same. There's room for optimization.

jchodera commented 7 years ago

OK, looks like I'm good!

I'll check this Abl:imatinib example into yank-examples when I've cleaned things up!

jchodera commented 7 years ago

https://github.com/choderalab/yank/pull/555 fixed some lingering issues!