I'm testing the PyAutoFEP and would like to start with the simplest case of morphing between two ligands.
I tried to use the star mapping scheme.

python ~/GitHub/PyAutoFEP/generate_perturbation_map.py --map_type=star --map_bias=RVXOH --input lig_data/*.mol2

However, I'm getting:

[ERROR] A total of molecules 2 were read. I need at least 3 to construct a meaningful perturbation graph
=================== STACK INFO ===================
  File "/Users/zhiyiwu/GitHub/PyAutoFEP/generate_perturbation_map.py", line 475, in <module>
    current_verbosity=arguments.verbose)
  File "/Users/zhiyiwu/GitHub/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================

I wonder if there is support for two molecules morphing as well?

If automated system builder is used, prepare_perturbation_map.py will fail if there are non-protein atoms in the input structure. Changes to prepare_perturbation_map.prepare_complex_system will be needed. The PDB order should be something like protein, cofactors/DNA/RNA, ligand, water, ions.

I´m stuck in the "Farnesoid X receptor tutorial" step in which the perturbation map is generated. I keep getting the following error:

Traceback (most recent call last):
File ".\PyAutoFEP\generate_perturbation_map.py", line 433, in
each_mol = mol_util.generic_mol_read(ligand_data=each_file, ligand_format=file_ext,
File ".\PyAutoFEP\mol_util.py", line 1211, in generic_mol_read
read_mol_data = rdkit.Chem.MolFromMolFile(ligand_data, removeHs=False)
OSError: Bad input file lig_data/*.mol

It seems that there's a problem with reading the files, so I also tried with simpler molecules and the mol2 file format, unfortunately getting the same "Bad input file" error.

I´m running the following command in a Conda environment:

python .\PyAutoFEP\generate_perturbation_map.py --map_type=star --map_bias=FXR_12 --input lig_data/*.mol

The Conda environment has installed the following packages:

python 3.9.12
rdkit 2022.03.3
networkx 2.8
alchemlyb 0.6.0
pymbar 3.0.5
openbabel 3.1.1
matplotlib 3.5.1
numpy 1.22.3

If anybody has any idea of how to fix this error, please let me know

SwissParam generates ligand topologies containing states A and B parameters. Currently, all_classes.TopologyData cannot parse such topologies.

• Add new __*data_fields and __data_dict to all_classes.TopologyData allowing for states A and B parameters.
• Add a fallback to each of add_bond, add_angle, add_pair, add_dihedral, and add_constraint, trying to read states A and B, checking if they are identical. NOTE: GROMACS manual [1] states that also position_restraints, dihedral_restraints, angle_restraints, and angle_restraints_z should accept states A and B, but none of this directives are processed by PyAutoFEP currently.
• Add a test using some SwissParam topologies

After the "Run FEP" step, the file "tutorial.tgz" has been created. However in the "nohup.out" file, the content are as follows:

tar: */water/md/rerun/*.pkl: Cannot stat: No such file or directory
tar: */water/md/rerun/*.xvg: Cannot stat: No such file or directory
tar: */water/md/rerun/*.edr: Cannot stat: No such file or directory
tar: */protein/md/rerun/*.pkl: Cannot stat: No such file or directory
tar: */protein/md/rerun/*.xvg: Cannot stat: No such file or directory
tar: */protein/md/rerun/*.edr: Cannot stat: No such file or directory
tar: */water/md/lambda*/lambda_pbcfit.xtc: Cannot stat: No such file or directory
tar: */water/md/lambda*/lambda.log: Cannot stat: No such file or directory
tar: */water/md/lambda*/lambda.edr: Cannot stat: No such file or directory
tar: */protein/md/lambda*/lambda_pbcfit.xtc: Cannot stat: No such file or directory
tar: */protein/md/lambda*/lambda.log: Cannot stat: No such file or directory
tar: */protein/md/lambda*/lambda.edr: Cannot stat: No such file or directory
tar: */water/md/lambda*/analysis.log: Cannot stat: No such file or directory
tar: */protein/md/lambda*/analysis.log: Cannot stat: No such file or directory
tar: */water/md/analysis/*.xvg: Cannot stat: No such file or directory
tar: */protein/md/analysis/*.xvg: Cannot stat: No such file or directory
tar: Exiting with failure status due to previous errors

And when i run the "analysis step", error happens,

[ERROR] Could not find unk_matrix.pkl file in tutorial/FXR_12-FXR_74/water/md/rerun to read temperature from. Also, could not s temperature from MD log in tutorial/FXR_12-FXR_74/water/md/rerun/../lambda0/lambda.log. Set temperature argument or config on. Analysis for this system will not be run.
=================== STACK INFO ===================
  File "analyze_results.py", line 1511, in <module>
    verbosity=arguments.verbose)
  File "analyze_results.py", line 1173, in collect_from_xvg
    msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
  File "/data/user/yanx/workdir/repo/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================

I think the 'rerun error' happens, there are no files in any 'rerun' fold. How should i do to fix this?
We would appreciate it if you could check with this when you have time.
@luancarvalhomartins

Edit: Formatted code output to improved readability - @luancarvalhomartins 2022.06.10

Dear Prof.

We set the following parameters:
"pose_loader = generic
poses_reference_structure = protein.pdb
prepare_dual_topology.py --poses_input $generic_poses --poses_input $generic_poses --config_file=input.ini"

The command works fine, but the results of pairs' overlap as if generic had not been added. It seems like "pose_loader = generic" isn't working.

We will appreciate your kindly advise.

Thanks,
Bao

================ Working on pairs ================
Perturbation Pose Coordinates
========== Working on ligand_1-ligand_2 ==========
[ERROR] Embeding of molecule ligand_2 (SMILES = [H]c1c(F)c(F)c(F)c([H])c1C([H])([H])n1c(=O)n(C([H])([H])C(=O)N([H])C([H])([H])[H])c(=O)n([H])/c1=N\c1c(C([H])([H])[H])c([H])c2nc(C([H])([H])[H])sc2c1[H]) to common core [10*]~c1~c~c~c~c~c~1~N~c1~n~c(~O)~n~c(~O)~n~1~C~c1~c~c(~F)~c(~F)~c(~F)~c~1~[11*] with molecule ligand_1 (SMILES = [H]c1c([H])c(OC([H])([H])[H])c([H])c([H])c1/N=c1\n([H])c(=O)n(C([H])([H])c2nc([H])n(C([H])([H])[H])n2)c(=O)n1C([H])([H])c1c([H])c(F)c(F)c(F)c1[H]) has failed. Molecule does not match the core.

We set up mcs_custom_mcs, and we made sure that SMARTS was the common core between Pairs. But we always get an error. We set the structure of the core to be a little bit larger.

We will appreciate your kindly advise.
Thanks,
Bao.

2022.05.27: I formatted the output for better reading (@luancarvalhomartins)

The following problem occurred while running the tutorial, and the previous steps were fine.

I am confused by the KeyError: 37, can not find any information of it.

(pyfep) jerome@ubuntu:~/Project/YXT/FEP-2/PyAutoFEP$ python3 generate_perturbation_map.py --map_type=star --map_bias=FXR_12 --input lig_data/*.mol
(pyfep) jerome@ubuntu:~/Project/YXT/FEP-2/PyAutoFEP$ python3 prepare_dual_topology.py --config_file=step2.ini
=============== mdp and run steps=================
         Complex                    Water          
        min01.mdp                 min01.mdp        
        min02.mdp                 min02.mdp        
        min03.mdp                  nve.mdp         
         nve.mdp                   nvt.mdp         
         nvt.mdp                   npt.mdp         
         npt.mdp                   md.mdp          
         md.mdp                       -            
==================================================
================= Input ligands ==================
   Name         Molecule            Topology      
  FXR_12   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Cl)c1[H] lig_data/FXR_12.itp
  FXR_74   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Br)c1[H] lig_data/FXR_74.itp
  FXR_76   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c([H])c1[H] lig_data/FXR_76.itp
  FXR_84   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(F)c1[H] lig_data/FXR_84.itp
  FXR_85   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C([H])([H])[H])c1[H] lig_data/FXR_85.itp
  FXR_88   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C(F)(F)F)c1[H] lig_data/FXR_88.itp
==================================================
=================== Poses read ===================
Name            File                
FXR_12          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d41c0>, 'topology': ['lig_data/FXR_12.itp']}
FXR_74          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4170>, 'topology': ['lig_data/FXR_74.itp']}
FXR_76          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4120>, 'topology': ['lig_data/FXR_76.itp']}
FXR_84          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4030>, 'topology': ['lig_data/FXR_84.itp']}
FXR_85          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfcd34e90>, 'topology': ['lig_data/FXR_85.itp']}
FXR_88          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4210>, 'topology': ['lig_data/FXR_88.itp']}
=============== Superimposed poses ===============
Name            File                      Note           
FXR_12          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d41c0>, 'topology': ['lig_data/FXR_12.itp']} Read from saved state
FXR_74          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4170>, 'topology': ['lig_data/FXR_74.itp']} Read from saved state
FXR_76          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4120>, 'topology': ['lig_data/FXR_76.itp']} Read from saved state
FXR_84          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4030>, 'topology': ['lig_data/FXR_84.itp']} Read from saved state
FXR_85          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfcd34e90>, 'topology': ['lig_data/FXR_85.itp']} Read from saved state
FXR_88          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7efcfc0d4210>, 'topology': ['lig_data/FXR_88.itp']} Read from saved state
==================================================
================== Align poses ===================
Molecule FXR_12 aligned
Molecule FXR_74 aligned
Molecule FXR_76 aligned
Molecule FXR_84 aligned
Molecule FXR_85 aligned
Molecule FXR_88 aligned
================= Perturbations ==================
         State A                  State B         
         FXR_12                   FXR_74          
         FXR_12                   FXR_76          
         FXR_12                   FXR_84          
         FXR_12                   FXR_85          
         FXR_12                   FXR_88          
==================================================
================ Working on pairs ================
Perturbation Pose Coordinates
============ Working on FXR_12-FXR_74 ============
Traceback (most recent call last):
  File "prepare_dual_topology.py", line 3451, in <module>
    verbosity=arguments.verbose)
  File "/home/jerome/Project/YXT/FEP-2/PyAutoFEP/merge_topologies.py", line 728, in merge_topologies
    atom_b_charge = topology2.molecules[0].atoms_dict[b_idx].q_e
KeyError: 37

2022.03.29 Formatted the original bug report @luancarvalhomartins

Hi, I was able to generate the perturbation map, however, the topology preparation is giving me error.
It says failed to read pose data. Please, I need your help in this regard. Thanks!

PyAutoFEP]$ python prepare_dual_topology.py --config_file=step2.ini
=============== mdp and run steps=================
         Complex                    Water          
        min01.mdp                 min01.mdp        
        min02.mdp                 min02.mdp        
        min03.mdp                  nve.mdp         
         nve.mdp                   nvt.mdp         
         nvt.mdp                   npt.mdp         
         npt.mdp                   md.mdp          
         md.mdp                       -            
==================================================
================= Input ligands ==================
   Name         Molecule            Topology      
  FXR_12   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Cl)c1[H] docs/tutorial01/workdir/lig_data/FXR_12.itp
  FXR_74   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Br)c1[H] docs/tutorial01/workdir/lig_data/FXR_74.itp
  FXR_76   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c([H])c1[H] docs/tutorial01/workdir/lig_data/FXR_76.itp
  FXR_84   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(F)c1[H] docs/tutorial01/workdir/lig_data/FXR_84.itp
  FXR_85   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C([H])([H])[H])c1[H] docs/tutorial01/workdir/lig_data/FXR_85.itp
  FXR_88   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C(F)(F)F)c1[H] docs/tutorial01/workdir/lig_data/FXR_88.itp
==================================================
=================== Poses read ===================
Name            File                
[ERROR] Failed to read pose data from docs/tutorial01/workdir/receptor_data/9mv.pdb
[ERROR] 
[ERROR] 
[ERROR] 
[ERROR] output_scripttype = slurm
[ERROR] output_collecttype = python with type .pdb
[ERROR] 
[ERROR] 
[ERROR] 
[ERROR] output_scripttype = slurm
[ERROR] output_collecttype = python
=================== STACK INFO ===================
  File "prepare_dual_topology.py", line 3250, in <module>
    verbosity=arguments.verbose, **arguments.poses_advanced_options)
  File "prepare_dual_topology.py", line 1324, in align_ligands
    save_state=save_state, verbosity=verbosity, **kwargs)
  File "/home/PyAutoFEP/docking_readers/superimpose_loader.py", line 80, in superimpose_poses
    verbosity=verbosity)['reference']
  File "/home/PyAutoFEP/docking_readers/generic_loader.py", line 60, in extract_docking_poses
    docking_mol_rd = generic_mol_read(ligand_format, each_mol['molecule'], verbosity=verbosity)
  File "/home/PyAutoFEP/docking_readers/generic_loader.py", line 209, in generic_mol_read
    msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
  File "/home/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================

(I formatted the output as inline code, this greatly improves readability - @luancarvalhomartins)

(PyAutoFEP) ldd@ldd:~/Desktop $ analyze_results.py --input tutorial.tgz --units kcal --output_uncompress_directory tutorial --center_molecule=FXR_12

All available analysis will be run
[ERROR] Could not find unk_matrix.pkl file in tutorial/FXR_12-FXR_76/protein/md/rerun to read temperature from. Also, could not guess temperature from MD log in tutorial/FXR_12-FXR_76/protein/md/rerun/../lambda0/lambda.log. Set temperature argument or config option. Analysis for this system will not be run.
=================== STACK INFO ===================
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1508, in <module>
    verbosity=arguments.verbose)
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1170, in collect_from_xvg
    msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
  File "/home/ldd/app/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================
[ERROR] Could not find unk_matrix.pkl file in tutorial/FXR_12-FXR_76/water/md/rerun to read temperature from. Also, could not guess temperature from MD log in tutorial/FXR_12-FXR_76/water/md/rerun/../lambda0/lambda.log. Set temperature argument or config option. Analysis for this system will not be run.
=================== STACK INFO ===================
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1508, in <module>
    verbosity=arguments.verbose)
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1170, in collect_from_xvg
    msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
  File "/home/ldd/app/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================
[ERROR] Could not find unk_matrix.pkl file in tutorial/FXR_12-FXR_74/protein/md/rerun to read temperature from. Also, could not guess temperature from MD log in tutorial/FXR_12-FXR_74/protein/md/rerun/../lambda0/lambda.log. Set temperature argument or config option. Analysis for this system will not be run.
=================== STACK INFO ===================
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1508, in <module>
    verbosity=arguments.verbose)
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1170, in collect_from_xvg
    msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
  File "/home/ldd/app/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================
[ERROR] Could not find unk_matrix.pkl file in tutorial/FXR_12-FXR_74/water/md/rerun to read temperature from. Also, could not guess temperature from MD log in tutorial/FXR_12-FXR_74/water/md/rerun/../lambda0/lambda.log. Set temperature argument or config option. Analysis for this system will not be run.
=================== STACK INFO ===================
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1508, in <module>
    verbosity=arguments.verbose)
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1170, in collect_from_xvg
    msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
  File "/home/ldd/app/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================
[ERROR] The leg systems saved on your progress file are ['protein', 'water'] while the systems read from the input file tutorial.tgz are []. Data mismatch, cannot continue.
=================== STACK INFO ===================
  File "/home/ldd/app/PyAutoFEP/analyze_results.py", line 1585, in <module>
    msg_verbosity=os_util.verbosity_level.error)
  File "/home/ldd/app/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================

We run analyze_results.py, got the above results. We will appreciate your kindly advise.

Formatted the output (@luancarvalhomartins 2022.05.04)

In some systems, running gmx grompp to prepare a .tpr to gmx genion in prepare_complex_system and prepare_water_system fails with:

[ERROR] Failed to run gmx grompp. Error code 1.
(.......)
[ERROR] Fatal error:
[ERROR] number of coordinates in coordinate file
[ERROR] (/tmp/XXXX/build_system_SSMMHH_DDMMYYYY/protein/systemsolvated_step4_SSMMHH_DDMMYYYY.gro, NNNNN) does not match topology
[ERROR] (/tmp/XXXX/protein/build_system_SSMMHH_DDMMYYYY/systemsolv_step4_SSMMHH_DDMMYYYY.
top, NNNN)

This seems to happen when editing of systemsolv_step4_SSMMHH_DDMMYYYY.
top by gmx solvate is somehow delayed or buffer in the /tmp filesystem.

mol_util.obmol_to_rwmol fails to convert back OBMol to RWMol. align_ligands should avoid using openbabel's rotate and translate. Do that either directly or using RDKit.

================== Align poses ===================
Traceback (most recent call last):
File "/home/ying/2t/PyAutoFEP-master/prepare_dual_topology.py", line 3250, in
verbosity=arguments.verbose, **arguments.poses_advanced_options)
File "/home/ying/2t/PyAutoFEP-master/prepare_dual_topology.py", line 1361, in align_ligands
ligand_mol = mol_util.obmol_to_rwmol(each_ligand, verbosity=verbosity)
File "/home/ying/2t/PyAutoFEP-master/mol_util.py", line 205, in obmol_to_rwmol
ignoreSmoothingFailures=True)
RuntimeError: Pre-condition Violation
getNumImplicitHs() called without preceding call to calcImplicitValence()
Violation occurred on line 183 in file Code/GraphMol/Atom.cpp
Failed Expression: d_implicitValence > -1
RDKIT: 2021.09.5
BOOST: 1_74

Discussed in #18

While I am processing the step below followed by tutorial:
$prepare_dual_topology.py --config_file=step2.ini
I got this error message:
=============== Superimposed poses ===============
Name File Note
Traceback (most recent call last):
File "/home/ying/2t/PyAutoFEP-master/prepare_dual_topology.py", line 3245, in
poses_mol_data = align_ligands(receptor_structure_mol, poses_input,
File "/home/ying/2t/PyAutoFEP-master/prepare_dual_topology.py", line 1323, in align_ligands
docking_mol_local = superimpose_poses(superimpose_loader_ligands, reference_pose_superimpose,
File "/home/ying/2t/PyAutoFEP-master/docking_readers/superimpose_loader.py", line 124, in superimpose_poses
thismol = merge_topologies.constrained_embed_shapeselect(each_ligand_mol, rdkit_reference_pose,
File "/home/ying/2t/PyAutoFEP-master/merge_topologies.py", line 360, in constrained_embed_shapeselect
[molecule.RemoveConformer(conformer.GetId()) for conformer in molecule.GetConformers()
File "/home/ying/2t/PyAutoFEP-master/merge_topologies.py", line 360, in
[molecule.RemoveConformer(conformer.GetId()) for conformer in molecule.GetConformers()
RuntimeError: Sequence modified during iteration

I'm sure rdkit and openbabel has been installed, only one point is that I changed 'import pybel' to 'from openbabel import pybel' in all the python scripts, how can I fix it? May you help me, thanks.

Our input molecules including imidazole fragments,. It seems like rdkit can't kekulize mol at Superimosed poses step. Is there any way to fix it?

================= Input ligands ==================
Name Molecule Topology
A [H]c1c([H])c(C(F)(F)F)c([H])c([H])c1N([H])c1c([H])c([H])c(S(=O)(=O)N([H])C([H])([H])[H])c([H])c1-c1nc([H])n(C([H])([H])[H])c1[H] lig_data/A.itp
B [H]c1c([H])c(C(F)(F)F)c([H])c([H])c1Oc1c([H])c([H])c(S(=O)(=O)N([H])C([H])([H])[H])c([H])c1-c1nc([H])n(C([H])([H])[H])c1[H] lig_data/B.itp
C [H]c1nc(N([H])c2c([H])c([H])c(S(=O)(=O)N([H])C([H])([H])[H])c([H])c2-c2nc([H])n(C([H])([H])[H])c2[H])c([H])c([H])c1C(F)(F)F lig_data/C.itp
==================================================
=================== Poses read ===================
Name File Details
A {'molecule': <rdkit.Chem.PropertyMol.PropertyMol object at 0x7f31f741dbf0>, 'topology': ['lig_data/A.itp']}
B {'molecule': <rdkit.Chem.PropertyMol.PropertyMol object at 0x7f323d985df0>, 'topology': ['lig_data/B.itp']}
C {'molecule': <rdkit.Chem.PropertyMol.PropertyMol object at 0x7f323e5e53b0>, 'topology': ['lig_data/C.itp']}
=============== Superimposed poses ===============
Name File Note
Traceback (most recent call last):
File "../prepare_dual_topology.py", line 3451, in
verbosity=arguments.verbose, **arguments.poses_advanced_options)
File "../prepare_dual_topology.py", line 2878, in align_ligands
save_state=save_state, verbosity=verbosity, **kwargs)
File "/home/zorrosword/PyAutoFEP/docking_readers/superimpose_loader.py", line 127, in superimpose_poses
atom_map=this_atommap, **kwargs)
File "/home/zorrosword/PyAutoFEP/os_util.py", line 517, in wrap_trace
return f(*args, **kwargs)
File "/home/zorrosword/PyAutoFEP/merge_topologies.py", line 243, in constrained_embed_shapeselect
core_mol = rdkit.Chem.RemoveHs(core_mol)
rdkit.Chem.rdchem.KekulizeException: Can't kekulize mol. Unkekulized atoms: 20

we will appreciate your kindly advise.

Thanks,

prepare_dual_topology.py --config_file=step2.ini --output_hidden_temp_dir=false
=============== mdp and run steps=================
         Complex                    Water          
        min01.mdp                 min01.mdp        
        min02.mdp                 min02.mdp        
        min03.mdp                  nve.mdp         
         nve.mdp                   nvt.mdp         
         nvt.mdp                   npt.mdp         
         npt.mdp                   md.mdp          
         md.mdp                       -            
==================================================
================= Input ligands ==================
   Name         Molecule            Topology      
  drg_1    [H]C#Cc1c(F)c([H])c([H])c2c([H])c(O[H])c([H])c(-c3nc([H])c4c(N5C([H])([H])C6([H])C([H])([H])C([H])([H])C([H])(C5([H])[H])[N+]6([H])[H])nc(OC([H])([H])C56C([H])([H])C([H])([H])C([H])([H])[N+]5([H])C([H])([H])C([H])(F)C6([H])[H])nc4c3Cl)c12 lig_data/drg_1.itp 
  drg_2    [H]C#Cc1c(Cl)c([H])c([H])c2c([H])c(O[H])c([H])c(-c3nc([H])c4c(N5C([H])([H])C6([H])C([H])([H])C([H])([H])C([H])(C5([H])[H])[N+]6([H])[H])nc(OC([H])([H])C56C([H])([H])C([H])([H])C([H])([H])[N+]5([H])C([H])([H])C([H])(F)C6([H])[H])nc4c3F)c12 lig_data/drg_2.itp 
  drg_3    [H]C#Cc1c(F)c([H])c([H])c2c([H])c([H])c([H])c(-c3nc([H])c4c(N5C([H])([H])C6([H])C([H])([H])C([H])([H])C([H])(C5([H])[H])[N+]6([H])[H])nc(OC([H])([H])C56C([H])([H])C([H])([H])C([H])([H])[N+]5([H])C([H])([H])C([H])(F)C6([H])[H])nc4c3F)c12 lig_data/drg_3.itp 
==================================================
=================== Poses read ===================
Name            File                 Details             
drg_1           {'molecule': <rdkit.Chem.PropertyMol.PropertyMol object at 0x7f52b40520d8>, 'topology': ['lig_data/drg_1.itp']}                   
drg_2           {'molecule': <rdkit.Chem.PropertyMol.PropertyMol object at 0x7f52b4052298>, 'topology': ['lig_data/drg_2.itp']}                   
drg_3           {'molecule': <rdkit.Chem.PropertyMol.PropertyMol object at 0x7f52b4052110>, 'topology': ['lig_data/drg_3.itp']}                   
=============== Superimposed poses ===============
Name            File                      Note           
Traceback (most recent call last):
  File "/home/cadd/program/fep/prepare_dual_topology.py", line 3451, in <module>
    verbosity=arguments.verbose, **arguments.poses_advanced_options)
  File "/home/cadd/program/fep/prepare_dual_topology.py", line 2878, in align_ligands
    save_state=save_state, verbosity=verbosity, **kwargs)
  File "/home/cadd/program/fep/docking_readers/superimpose_loader.py", line 127, in superimpose_poses
    atom_map=this_atommap, **kwargs)
  File "/home/cadd/program/fep/os_util.py", line 517, in wrap_trace
    return f(*args, **kwargs)
  File "/home/cadd/program/fep/merge_topologies.py", line 243, in constrained_embed_shapeselect
    core_mol = rdkit.Chem.RemoveHs(core_mol)
rdkit.Chem.rdchem.MolSanitizeException: Explicit valence for atom # 7 C, 5, is greater than permitted

workdir.zip

2022-07-12: formatting (@luancarvalhomartins)

Gromacs2022.1 has been released and I ran the sample and got the following error. There is nothing wrong with the file or other files. When I used the same files and ran them with Gromacs2020.6, it worked fine. We would appreciate it if you could check with 2022.1 when you have time.

[ERROR] Failed to run gmx pdb2gmx. Error code 1.
[ERROR] Command line was: ['gmx', 'pdb2gmx', '-f', 'receptor_data/5q17_processed.pdb', '-p', '/work113/tmp/PyAutoFEP/docs/tutorial01_plumed2022/tutorial/FXR_12-FXR_74/protein/build_system_125508_26042022/protein_step1_125508_26042022.top', '-o', '/work113/tmp/PyAutoFEP/docs/tutorial01_plumed2022/tutorial/FXR_12-FXR_74/protein/build_system_125508_26042022/protein_step1_125508_26042022.pdb']
[ERROR]
[ERROR] stdout:
[ERROR] Select the Force Field:
[ERROR]
[ERROR] From current directory:
[ERROR]
[ERROR]  1: OPLS-AA/M all-atom force field (2015 aminoacid dihedrals)
[ERROR]
[ERROR] From '/home/gromacs-2022.1_plumed/share/gromacs/top':
[ERROR]
[ERROR]  2: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
[ERROR]
[ERROR]  3: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
[ERROR]
[ERROR]  4: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
[ERROR]
[ERROR]  5: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
[ERROR]
[ERROR]  6: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
[ERROR]
[ERROR]  7: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
[ERROR]
[ERROR]  8: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
[ERROR]
[ERROR]  9: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
[ERROR]
[ERROR] 10: GROMOS96 43a1 force field
[ERROR]
[ERROR] 11: GROMOS96 43a2 force field (improved alkane dihedrals)
[ERROR]
[ERROR] 12: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
[ERROR]
[ERROR] 13: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
[ERROR]
[ERROR] 14: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
[ERROR]
[ERROR] 15: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
[ERROR]
[ERROR] 16: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
[ERROR]
[ERROR] Using the Oplsaam force field in directory ./oplsaam.ff
[ERROR]
[ERROR] Select the Water Model:
[ERROR]
[ERROR]  1: TIP4P  TIP 4-point, recommended
[ERROR]
[ERROR]  2: TIP3P  TIP 3-point
[ERROR]
[ERROR]  3: TIP5P  TIP 5-point (see http://redmine.gromacs.org/issues/1348 for issues)
[ERROR]
[ERROR]  4: TIP5P  TIP 5-point improved for Ewald sums
[ERROR]
[ERROR]  5: SPC    simple point charge
[ERROR]
[ERROR]  6: SPC/E  extended simple point charge
[ERROR]
[ERROR]  7: None
[ERROR]
[ERROR] going to rename ./oplsaam.ff/aminoacids.r2b
[ERROR] Reading receptor_data/5q17_processed.pdb...
[ERROR] Read '', 3871 atoms
[ERROR]
[ERROR] Analyzing pdb file
[ERROR] Splitting chemical chains based on TER records or chain id changing.
[ERROR]
[ERROR] There are 2 chains and 0 blocks of water and 236 residues with 3871 atoms
[ERROR]
[ERROR]   chain  #res #atoms
[ERROR]
[ERROR]   1 ' '   224   3662
[ERROR]
[ERROR]   2 ' '    12    209
[ERROR]
[ERROR] there were 2106 atoms with zero occupancy and 1765 atoms with          occupancy unequal to one (out of 3871 atoms). Check your pdb file.
[ERROR]
[ERROR] Reading residue database... (Oplsaam)
[ERROR]
[ERROR]
[ERROR] stderr:
[ERROR]                :-) GROMACS - gmx pdb2gmx, 2022.1-plumed_2.8.0 (-:
[ERROR]
[ERROR] Executable:   /home/gromacs-2022.1_plumed/bin/gmx
[ERROR] Data prefix:  /home/gromacs-2022.1_plumed
[ERROR] Working dir:  /work113/tmp/PyAutoFEP/docs/tutorial01_plumed2022
[ERROR] Command line:
[ERROR]   gmx pdb2gmx -f receptor_data/5q17_processed.pdb -p /work113/tmp/PyAutoFEP/docs/tutorial01_plumed2022/tutorial/FXR_12-FXR_74/protein/build_system_125508_26042022/protein_step1_125508_26042022.top -o /work113/tmp/PyAutoFEP/docs/tutorial01_plumed2022/tutorial/FXR_12-FXR_74/protein/build_system_125508_26042022/protein_step1_125508_26042022.pdb
[ERROR]
[ERROR] Opening force field file ./oplsaam.ff/watermodels.dat
[ERROR] Opening force field file ./oplsaam.ff/aminoacids.r2b
[ERROR] there were 2106 atoms with zero occupancy and 1765 atoms with          occupancy unequal to one (out of 3871 atoms). Check your pdb file.
[ERROR] Opening force field file ./oplsaam.ff/atomtypes.atp
[ERROR] Opening force field file ./oplsaam.ff/aminoacids.rtp
[ERROR] Opening force field file ./oplsaam.ff/aminoacids.hdb
[ERROR] Opening force field file ./oplsaam.ff/aminoacids.n.tdb
[ERROR]
[ERROR] -------------------------------------------------------
[ERROR] Program:     gmx pdb2gmx, version 2022.1-plumed_2.8.0
[ERROR] Source file: src/gromacs/gmxpreprocess/ter_db.cpp (line 139)
[ERROR] Function:    void read_atom(char*, bool, std::__cxx11::string*, t_atom*, PreprocessingAtomTypes*, int*)
[ERROR]
[ERROR] Inconsistency in user input:
[ERROR] Atomtype for atom name oplsm_298 not found in terminal data base
[ERROR]
[ERROR] For more information and tips for troubleshooting, please check the GROMACS
[ERROR] website at http://www.gromacs.org/Documentation/Errors
[ERROR] -------------------------------------------------------
[ERROR]
=================== STACK INFO ===================
  File "../../prepare_dual_topology.py", line 3445, in <module>
    solvate_data=solvate_data, verbosity=arguments.verbose)
  File "../../prepare_dual_topology.py", line 249, in prepare_complex_system
    os_util.run_gmx(gmx_bin, pdb2gmx_list, communicate_str, build_files_dict['pdb2gmx_log'], verbosity=verbosity)
  File "/work113/tmp/PyAutoFEP/os_util.py", line 146, in run_gmx
    msg_verbosity=verbosity_level.error, current_verbosity=verbosity)
  File "/work113/tmp/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================

Originally posted by @baba-hashimoto in #30 (comment)

Thanks for sharing your nice works, could you create a jupyter notebook to describe how to use your scripts? It will be a great help for new guys to PyAutoFEP

I set "extradirs = /usr/local/gromacs/share/gromacs/top/amber99sb-ildn.ff", it doesn't work.

When i run the tutorial step by step and got the following error. We would appreciate it if you could check with the tutorial with Gromacs 2020 when you have time.

[ERROR] Failed to run gmx grompp. Error code 1.
[ERROR] Command line was: ['gmx', 'grompp', '-f', '/tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/genion_step5_160215_08062022.mdp', '-c', '/tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/systemsolvated_step4_160215_08062022.gro', '-p', '/tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/systemsolv_step4_160215_08062022.top', '-o', '/tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/genion_step6_160215_08062022.tpr', '-maxwarn', '1', '-po', '/tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/mdout_step6_160215_08062022.mdp']
[ERROR]
[ERROR] stdout:
[ERROR]
[ERROR]
[ERROR] stderr:
[ERROR] :-) GROMACS - gmx grompp, 2020.6 (-:
[ERROR]
[ERROR] GROMACS is written by:
[ERROR] Emile Apol Rossen Apostolov Paul Bauer Herman J.C. Berendsen
[ERROR] Par Bjelkmar Christian Blau Viacheslav Bolnykh Kevin Boyd
[ERROR] Aldert van Buuren Rudi van Drunen Anton Feenstra Alan Gray
[ERROR] Gerrit Groenhof Anca Hamuraru Vincent Hindriksen M. Eric Irrgang
[ERROR] Aleksei Iupinov Christoph Junghans Joe Jordan Dimitrios Karkoulis
[ERROR] Peter Kasson Jiri Kraus Carsten Kutzner Per Larsson
[ERROR] Justin A. Lemkul Viveca Lindahl Magnus Lundborg Erik Marklund
[ERROR] Pascal Merz Pieter Meulenhoff Teemu Murtola Szilard Pall
[ERROR] Sander Pronk Roland Schulz Michael Shirts Alexey Shvetsov
[ERROR] Alfons Sijbers Peter Tieleman Jon Vincent Teemu Virolainen
[ERROR] Christian Wennberg Maarten Wolf Artem Zhmurov
[ERROR] and the project leaders:
[ERROR] Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel
[ERROR]
[ERROR] Copyright (c) 1991-2000, University of Groningen, The Netherlands.
[ERROR] Copyright (c) 2001-2019, The GROMACS development team at
[ERROR] Uppsala University, Stockholm University and
[ERROR] the Royal Institute of Technology, Sweden.
[ERROR] check out http://www.gromacs.org for more information.
[ERROR]
[ERROR] GROMACS is free software; you can redistribute it and/or modify it
[ERROR] under the terms of the GNU Lesser General Public License
[ERROR] as published by the Free Software Foundation; either version 2.1
[ERROR] of the License, or (at your option) any later version.
[ERROR]
[ERROR] GROMACS: gmx grompp, version 2020.6
[ERROR] Executable: /data/software/gmx2020.6/bin/gmx
[ERROR] Data prefix: /data/software/gmx2020.6
[ERROR] Working dir: /data/user/yanx/workdir/repo/PyAutoFEP
[ERROR] Command line:
[ERROR] gmx grompp -f /tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/genion_step5_160215_08062022.mdp -c /tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/systemsolvated_step4_160215_08062022.gro -p /tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/systemsolv_step4_160215_08062022.top -o /tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/genion_step6_160215_08062022.tpr -maxwarn 1 -po /tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/mdout_step6_160215_08062022.mdp
[ERROR]
[ERROR]
[ERROR] NOTE 1 [file /tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/genion_step5_160215_08062022.mdp]:
[ERROR] For a correct single-point energy evaluation with nsteps = 0, use
[ERROR] continuation = yes to avoid constraining the input coordinates.
[ERROR]
[ERROR] Setting the LD random seed to -44077602
[ERROR] Generated 449826 of the 449826 non-bonded parameter combinations
[ERROR] Generating 1-4 interactions: fudge = 0.5
[ERROR] Generated 449826 of the 449826 1-4 parameter combinations
[ERROR] Excluding 3 bonded neighbours molecule type 'Protein'
[ERROR] Excluding 3 bonded neighbours molecule type 'Protein2'
[ERROR] Excluding 3 bonded neighbours molecule type 'LIG'
[ERROR]
[ERROR] NOTE 2 [file systemsolv_step4_160215_08062022.top, line 56]:
[ERROR] System has non-zero total charge: -10.000000
[ERROR] Total charge should normally be an integer. See
[ERROR] http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
[ERROR] for discussion on how close it should be to an integer.
[ERROR]
[ERROR]
[ERROR]
[ERROR]
[ERROR] -------------------------------------------------------
[ERROR] Program: gmx grompp, version 2020.6
[ERROR] Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 669)
[ERROR]
[ERROR] Fatal error:
[ERROR] number of coordinates in coordinate file
[ERROR] (/tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/systemsolvated_step4_160215_08062022.gro,
[ERROR] 44742)
[ERROR] does not match topology
[ERROR] (/tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/systemsolv_step4_160215_08062022.top,
[ERROR] 3930)
[ERROR]
[ERROR] For more information and tips for troubleshooting, please check the GROMACS
[ERROR] website at http://www.gromacs.org/Documentation/Errors
[ERROR] -------------------------------------------------------
[ERROR]
[ERROR]
[ERROR] ==================================================
[ERROR] This is likely caused by a failing to edit the intermediate topology file /tmp/tmpmxflrejy/tutorial/FXR_12-FXR_74/protein/build_system_160215_08062022/systemsolv_step4_160215_08062022.top. Rerunning with output_hidden_temp_dir=False may solve this issue.
=================== STACK INFO ===================
File "prepare_dual_topology.py", line 3547, in
verbosity=arguments.verbose)
File "prepare_dual_topology.py", line 494, in prepare_complex_system
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
File "/data/user/yanx/workdir/repo/PyAutoFEP/os_util.py", line 292, in local_print
formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO =================== @luancarvalhomartins

I'm interested in understanding how the topology files are made.
I have gone through the tutorial and the tutorial.bin would extract to a number of legs.
For example, for FXR_12-FXR_74/protein/lambda4/ligand.itp.

It looks like

[ moleculetype ]
; Name               nrexcl
LIG        3
[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  
1       O00_const_idx1 1       LIG     O00_const_idx1 1       -0.607809090909091 15.999  
2       C01_const_idx2 1       LIG     C01_const_idx2 1       0.4391090909090909 12.011  
3       O02_const_idx3 1       LIG     O02_const_idx3 1       -0.607809090909091 15.999  

I was kind of expecting a file that looks likes https://manual.gromacs.org/documentation/current/reference-manual/topologies/topology-file-formats.html#topologies-for-free-energy-calculations

[ moleculetype ]
; Name               nrexcl
LIG        3
[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  
1       O00_const_idx1 1       LIG     O00_const_idx1 1       -0.607809090909091 15.999  O00_const_idx1 -0 15.999
2       C01_const_idx2 1       LIG     C01_const_idx2 1       0.4391090909090909 12.011  C01_const_idx2 -0 12.011
3       O02_const_idx3 1       LIG     O02_const_idx3 1       -0.607809090909091 15.999  O02_const_idx3 -0 15.999

Where B state is the end state. I wonder if you mind helping me understand the file setup, please? Thank you.

Hi,
I'm interested in using PyAutoFEP for my project. For reproducibility, I wonder if there is a stable release that is recommended?
Thank you.

Hi. I still encountered environment: rerun/rerun_struct_coord.out: No such file or directory errors.
My perturbation_dir is /program/fep/workdir/fep, and I get the unk_matrix.pkl file in rerun/
unk_matrix.zip

cadd@huiyu:~/program/fep/workdir/fep/drg_1-drg_2/protein/md/rerun$ ls

lambda0.xvg rerun_struct11_coord8.edr rerun_struct3_coord9.edr rerun_struct7_coord0.edr lambda10.xvg rerun_struct11_coord8.log rerun_struct3_coord9.log rerun_struct7_coord0.log lambda11.xvg rerun_struct11_coord8.out rerun_struct3_coord9.out rerun_struct7_coord0.out lambda1.xvg rerun_struct11_coord8.xvg rerun_struct3_coord9.xvg rerun_struct7_coord0.xvg lambda2.xvg rerun_struct11_coord9.edr rerun_struct4_coord0.edr rerun_struct7_coord10.edr lambda3.xvg rerun_struct11_coord9.log rerun_struct4_coord0.log rerun_struct7_coord10.log lambda4.xvg rerun_struct11_coord9.out rerun_struct4_coord0.out rerun_struct7_coord10.out lambda5.xvg rerun_struct11_coord9.xvg rerun_struct4_coord0.xvg rerun_struct7_coord10.xvg lambda6.xvg rerun_struct1_coord0.edr rerun_struct4_coord10.edr rerun_struct7_coord11.edr lambda7.xvg rerun_struct1_coord0.log rerun_struct4_coord10.log rerun_struct7_coord11.log lambda8.xvg rerun_struct1_coord0.out rerun_struct4_coord10.out rerun_struct7_coord11.out lambda9.xvg rerun_struct1_coord0.xvg rerun_struct4_coord10.xvg rerun_struct7_coord11.xvg rerun_struct0_coord0.edr rerun_struct1_coord10.edr rerun_struct4_coord11.edr rerun_struct7_coord1.edr rerun_struct0_coord0.log rerun_struct1_coord10.log rerun_struct4_coord11.log rerun_struct7_coord1.log rerun_struct0_coord0.out rerun_struct1_coord10.out rerun_struct4_coord11.out rerun_struct7_coord1.out rerun_struct0_coord0.xvg rerun_struct1_coord10.xvg rerun_struct4_coord11.xvg rerun_struct7_coord1.xvg rerun_struct0_coord10.edr rerun_struct1_coord11.edr rerun_struct4_coord1.edr rerun_struct7_coord2.edr rerun_struct0_coord10.log rerun_struct1_coord11.log rerun_struct4_coord1.log rerun_struct7_coord2.log rerun_struct0_coord10.out rerun_struct1_coord11.out rerun_struct4_coord1.out rerun_struct7_coord2.out rerun_struct0_coord10.xvg rerun_struct1_coord11.xvg rerun_struct4_coord1.xvg rerun_struct7_coord2.xvg rerun_struct0_coord11.edr rerun_struct1_coord1.edr rerun_struct4_coord2.edr rerun_struct7_coord3.edr rerun_struct0_coord11.log rerun_struct1_coord1.log rerun_struct4_coord2.log rerun_struct7_coord3.log rerun_struct0_coord11.out rerun_struct1_coord1.out rerun_struct4_coord2.out rerun_struct7_coord3.out rerun_struct0_coord11.xvg rerun_struct1_coord1.xvg rerun_struct4_coord2.xvg rerun_struct7_coord3.xvg rerun_struct0_coord1.edr rerun_struct1_coord2.edr rerun_struct4_coord3.edr rerun_struct7_coord4.edr rerun_struct0_coord1.log rerun_struct1_coord2.log rerun_struct4_coord3.log rerun_struct7_coord4.log rerun_struct0_coord1.out rerun_struct1_coord2.out rerun_struct4_coord3.out rerun_struct7_coord4.out rerun_struct0_coord1.xvg rerun_struct1_coord2.xvg rerun_struct4_coord3.xvg rerun_struct7_coord4.xvg rerun_struct0_coord2.edr rerun_struct1_coord3.edr rerun_struct4_coord4.edr rerun_struct7_coord5.edr rerun_struct0_coord2.log rerun_struct1_coord3.log rerun_struct4_coord4.log rerun_struct7_coord5.log rerun_struct0_coord2.out rerun_struct1_coord3.out rerun_struct4_coord4.out rerun_struct7_coord5.out rerun_struct0_coord2.xvg rerun_struct1_coord3.xvg rerun_struct4_coord4.xvg rerun_struct7_coord5.xvg rerun_struct0_coord3.edr rerun_struct1_coord4.edr rerun_struct4_coord5.edr rerun_struct7_coord6.edr rerun_struct0_coord3.log rerun_struct1_coord4.log rerun_struct4_coord5.log rerun_struct7_coord6.log rerun_struct0_coord3.out rerun_struct1_coord4.out rerun_struct4_coord5.out rerun_struct7_coord6.out rerun_struct0_coord3.xvg rerun_struct1_coord4.xvg rerun_struct4_coord5.xvg rerun_struct7_coord6.xvg rerun_struct0_coord4.edr rerun_struct1_coord5.edr rerun_struct4_coord6.edr rerun_struct7_coord7.edr rerun_struct0_coord4.log rerun_struct1_coord5.log rerun_struct4_coord6.log rerun_struct7_coord7.log rerun_struct0_coord4.out rerun_struct1_coord5.out rerun_struct4_coord6.out rerun_struct7_coord7.out rerun_struct0_coord4.xvg rerun_struct1_coord5.xvg rerun_struct4_coord6.xvg rerun_struct7_coord7.xvg rerun_struct0_coord5.edr rerun_struct1_coord6.edr rerun_struct4_coord7.edr rerun_struct7_coord8.edr rerun_struct0_coord5.log rerun_struct1_coord6.log rerun_struct4_coord7.log rerun_struct7_coord8.log rerun_struct0_coord5.out rerun_struct1_coord6.out rerun_struct4_coord7.out rerun_struct7_coord8.out rerun_struct0_coord5.xvg rerun_struct1_coord6.xvg rerun_struct4_coord7.xvg rerun_struct7_coord8.xvg rerun_struct0_coord6.edr rerun_struct1_coord7.edr rerun_struct4_coord8.edr rerun_struct7_coord9.edr rerun_struct0_coord6.log rerun_struct1_coord7.log rerun_struct4_coord8.log rerun_struct7_coord9.log rerun_struct0_coord6.out rerun_struct1_coord7.out rerun_struct4_coord8.out rerun_struct7_coord9.out rerun_struct0_coord6.xvg rerun_struct1_coord7.xvg rerun_struct4_coord8.xvg rerun_struct7_coord9.xvg rerun_struct0_coord7.edr rerun_struct1_coord8.edr rerun_struct4_coord9.edr rerun_struct8_coord0.edr rerun_struct0_coord7.log rerun_struct1_coord8.log rerun_struct4_coord9.log rerun_struct8_coord0.log rerun_struct0_coord7.out rerun_struct1_coord8.out rerun_struct4_coord9.out rerun_struct8_coord0.out rerun_struct0_coord7.xvg rerun_struct1_coord8.xvg rerun_struct4_coord9.xvg rerun_struct8_coord0.xvg rerun_struct0_coord8.edr rerun_struct1_coord9.edr rerun_struct5_coord0.edr rerun_struct8_coord10.edr rerun_struct0_coord8.log rerun_struct1_coord9.log rerun_struct5_coord0.log rerun_struct8_coord10.log rerun_struct0_coord8.out rerun_struct1_coord9.out rerun_struct5_coord0.out rerun_struct8_coord10.out rerun_struct0_coord8.xvg rerun_struct1_coord9.xvg rerun_struct5_coord0.xvg rerun_struct8_coord10.xvg rerun_struct0_coord9.edr rerun_struct2_coord0.edr rerun_struct5_coord10.edr rerun_struct8_coord11.edr rerun_struct0_coord9.log rerun_struct2_coord0.log rerun_struct5_coord10.log rerun_struct8_coord11.log rerun_struct0_coord9.out rerun_struct2_coord0.out rerun_struct5_coord10.out rerun_struct8_coord11.out rerun_struct0_coord9.xvg rerun_struct2_coord0.xvg rerun_struct5_coord10.xvg rerun_struct8_coord11.xvg rerun_struct10_coord0.edr rerun_struct2_coord10.edr rerun_struct5_coord11.edr rerun_struct8_coord1.edr rerun_struct10_coord0.log rerun_struct2_coord10.log rerun_struct5_coord11.log rerun_struct8_coord1.log rerun_struct10_coord0.out rerun_struct2_coord10.out rerun_struct5_coord11.out rerun_struct8_coord1.out rerun_struct10_coord0.xvg rerun_struct2_coord10.xvg rerun_struct5_coord11.xvg rerun_struct8_coord1.xvg rerun_struct10_coord10.edr rerun_struct2_coord11.edr rerun_struct5_coord1.edr rerun_struct8_coord2.edr rerun_struct10_coord10.log rerun_struct2_coord11.log rerun_struct5_coord1.log rerun_struct8_coord2.log rerun_struct10_coord10.out rerun_struct2_coord11.out rerun_struct5_coord1.out rerun_struct8_coord2.out rerun_struct10_coord10.xvg rerun_struct2_coord11.xvg rerun_struct5_coord1.xvg rerun_struct8_coord2.xvg rerun_struct10_coord11.edr rerun_struct2_coord1.edr rerun_struct5_coord2.edr rerun_struct8_coord3.edr rerun_struct10_coord11.log rerun_struct2_coord1.log rerun_struct5_coord2.log rerun_struct8_coord3.log rerun_struct10_coord11.out rerun_struct2_coord1.out rerun_struct5_coord2.out rerun_struct8_coord3.out rerun_struct10_coord11.xvg rerun_struct2_coord1.xvg rerun_struct5_coord2.xvg rerun_struct8_coord3.xvg rerun_struct10_coord1.edr rerun_struct2_coord2.edr rerun_struct5_coord3.edr rerun_struct8_coord4.edr rerun_struct10_coord1.log rerun_struct2_coord2.log rerun_struct5_coord3.log rerun_struct8_coord4.log rerun_struct10_coord1.out rerun_struct2_coord2.out rerun_struct5_coord3.out rerun_struct8_coord4.out rerun_struct10_coord1.xvg rerun_struct2_coord2.xvg rerun_struct5_coord3.xvg rerun_struct8_coord4.xvg rerun_struct10_coord2.edr rerun_struct2_coord3.edr rerun_struct5_coord4.edr rerun_struct8_coord5.edr rerun_struct10_coord2.log rerun_struct2_coord3.log rerun_struct5_coord4.log rerun_struct8_coord5.log rerun_struct10_coord2.out rerun_struct2_coord3.out rerun_struct5_coord4.out rerun_struct8_coord5.out rerun_struct10_coord2.xvg rerun_struct2_coord3.xvg rerun_struct5_coord4.xvg rerun_struct8_coord5.xvg rerun_struct10_coord3.edr rerun_struct2_coord4.edr rerun_struct5_coord5.edr rerun_struct8_coord6.edr rerun_struct10_coord3.log rerun_struct2_coord4.log rerun_struct5_coord5.log rerun_struct8_coord6.log rerun_struct10_coord3.out rerun_struct2_coord4.out rerun_struct5_coord5.out rerun_struct8_coord6.out rerun_struct10_coord3.xvg rerun_struct2_coord4.xvg rerun_struct5_coord5.xvg rerun_struct8_coord6.xvg rerun_struct10_coord4.edr rerun_struct2_coord5.edr rerun_struct5_coord6.edr rerun_struct8_coord7.edr rerun_struct10_coord4.log rerun_struct2_coord5.log rerun_struct5_coord6.log rerun_struct8_coord7.log rerun_struct10_coord4.out rerun_struct2_coord5.out rerun_struct5_coord6.out rerun_struct8_coord7.out rerun_struct10_coord4.xvg rerun_struct2_coord5.xvg rerun_struct5_coord6.xvg rerun_struct8_coord7.xvg rerun_struct10_coord5.edr rerun_struct2_coord6.edr rerun_struct5_coord7.edr rerun_struct8_coord8.edr rerun_struct10_coord5.log rerun_struct2_coord6.log rerun_struct5_coord7.log rerun_struct8_coord8.log rerun_struct10_coord5.out rerun_struct2_coord6.out rerun_struct5_coord7.out rerun_struct8_coord8.out rerun_struct10_coord5.xvg rerun_struct2_coord6.xvg rerun_struct5_coord7.xvg rerun_struct8_coord8.xvg rerun_struct10_coord6.edr rerun_struct2_coord7.edr rerun_struct5_coord8.edr rerun_struct8_coord9.edr rerun_struct10_coord6.log rerun_struct2_coord7.log rerun_struct5_coord8.log rerun_struct8_coord9.log rerun_struct10_coord6.out rerun_struct2_coord7.out rerun_struct5_coord8.out rerun_struct8_coord9.out rerun_struct10_coord6.xvg rerun_struct2_coord7.xvg rerun_struct5_coord8.xvg rerun_struct8_coord9.xvg rerun_struct10_coord7.edr rerun_struct2_coord8.edr rerun_struct5_coord9.edr rerun_struct9_coord0.edr rerun_struct10_coord7.log rerun_struct2_coord8.log rerun_struct5_coord9.log rerun_struct9_coord0.log rerun_struct10_coord7.out rerun_struct2_coord8.out rerun_struct5_coord9.out rerun_struct9_coord0.out rerun_struct10_coord7.xvg rerun_struct2_coord8.xvg rerun_struct5_coord9.xvg rerun_struct9_coord0.xvg rerun_struct10_coord8.edr rerun_struct2_coord9.edr rerun_struct6_coord0.edr rerun_struct9_coord10.edr rerun_struct10_coord8.log rerun_struct2_coord9.log rerun_struct6_coord0.log rerun_struct9_coord10.log rerun_struct10_coord8.out rerun_struct2_coord9.out rerun_struct6_coord0.out rerun_struct9_coord10.out rerun_struct10_coord8.xvg rerun_struct2_coord9.xvg rerun_struct6_coord0.xvg rerun_struct9_coord10.xvg rerun_struct10_coord9.edr rerun_struct3_coord0.edr rerun_struct6_coord10.edr rerun_struct9_coord11.edr rerun_struct10_coord9.log rerun_struct3_coord0.log rerun_struct6_coord10.log rerun_struct9_coord11.log rerun_struct10_coord9.out rerun_struct3_coord0.out rerun_struct6_coord10.out rerun_struct9_coord11.out rerun_struct10_coord9.xvg rerun_struct3_coord0.xvg rerun_struct6_coord10.xvg rerun_struct9_coord11.xvg rerun_struct11_coord0.edr rerun_struct3_coord10.edr rerun_struct6_coord11.edr rerun_struct9_coord1.edr rerun_struct11_coord0.log rerun_struct3_coord10.log rerun_struct6_coord11.log rerun_struct9_coord1.log rerun_struct11_coord0.out rerun_struct3_coord10.out rerun_struct6_coord11.out rerun_struct9_coord1.out rerun_struct11_coord0.xvg rerun_struct3_coord10.xvg rerun_struct6_coord11.xvg rerun_struct9_coord1.xvg rerun_struct11_coord10.edr rerun_struct3_coord11.edr rerun_struct6_coord1.edr rerun_struct9_coord2.edr rerun_struct11_coord10.log rerun_struct3_coord11.log rerun_struct6_coord1.log rerun_struct9_coord2.log rerun_struct11_coord10.out rerun_struct3_coord11.out rerun_struct6_coord1.out rerun_struct9_coord2.out rerun_struct11_coord10.xvg rerun_struct3_coord11.xvg rerun_struct6_coord1.xvg rerun_struct9_coord2.xvg rerun_struct11_coord11.edr rerun_struct3_coord1.edr rerun_struct6_coord2.edr rerun_struct9_coord3.edr rerun_struct11_coord11.log rerun_struct3_coord1.log rerun_struct6_coord2.log rerun_struct9_coord3.log rerun_struct11_coord11.out rerun_struct3_coord1.out rerun_struct6_coord2.out rerun_struct9_coord3.out rerun_struct11_coord11.xvg rerun_struct3_coord1.xvg rerun_struct6_coord2.xvg rerun_struct9_coord3.xvg rerun_struct11_coord1.edr rerun_struct3_coord2.edr rerun_struct6_coord3.edr rerun_struct9_coord4.edr rerun_struct11_coord1.log rerun_struct3_coord2.log rerun_struct6_coord3.log rerun_struct9_coord4.log rerun_struct11_coord1.out rerun_struct3_coord2.out rerun_struct6_coord3.out rerun_struct9_coord4.out rerun_struct11_coord1.xvg rerun_struct3_coord2.xvg rerun_struct6_coord3.xvg rerun_struct9_coord4.xvg rerun_struct11_coord2.edr rerun_struct3_coord3.edr rerun_struct6_coord4.edr rerun_struct9_coord5.edr rerun_struct11_coord2.log rerun_struct3_coord3.log rerun_struct6_coord4.log rerun_struct9_coord5.log rerun_struct11_coord2.out rerun_struct3_coord3.out rerun_struct6_coord4.out rerun_struct9_coord5.out rerun_struct11_coord2.xvg rerun_struct3_coord3.xvg rerun_struct6_coord4.xvg rerun_struct9_coord5.xvg rerun_struct11_coord3.edr rerun_struct3_coord4.edr rerun_struct6_coord5.edr rerun_struct9_coord6.edr rerun_struct11_coord3.log rerun_struct3_coord4.log rerun_struct6_coord5.log rerun_struct9_coord6.log rerun_struct11_coord3.out rerun_struct3_coord4.out rerun_struct6_coord5.out rerun_struct9_coord6.out rerun_struct11_coord3.xvg rerun_struct3_coord4.xvg rerun_struct6_coord5.xvg rerun_struct9_coord6.xvg rerun_struct11_coord4.edr rerun_struct3_coord5.edr rerun_struct6_coord6.edr rerun_struct9_coord7.edr rerun_struct11_coord4.log rerun_struct3_coord5.log rerun_struct6_coord6.log rerun_struct9_coord7.log rerun_struct11_coord4.out rerun_struct3_coord5.out rerun_struct6_coord6.out rerun_struct9_coord7.out rerun_struct11_coord4.xvg rerun_struct3_coord5.xvg rerun_struct6_coord6.xvg rerun_struct9_coord7.xvg rerun_struct11_coord5.edr rerun_struct3_coord6.edr rerun_struct6_coord7.edr rerun_struct9_coord8.edr rerun_struct11_coord5.log rerun_struct3_coord6.log rerun_struct6_coord7.log rerun_struct9_coord8.log rerun_struct11_coord5.out rerun_struct3_coord6.out rerun_struct6_coord7.out rerun_struct9_coord8.out rerun_struct11_coord5.xvg rerun_struct3_coord6.xvg rerun_struct6_coord7.xvg rerun_struct9_coord8.xvg rerun_struct11_coord6.edr rerun_struct3_coord7.edr rerun_struct6_coord8.edr rerun_struct9_coord9.edr rerun_struct11_coord6.log rerun_struct3_coord7.log rerun_struct6_coord8.log rerun_struct9_coord9.log rerun_struct11_coord6.out rerun_struct3_coord7.out rerun_struct6_coord8.out rerun_struct9_coord9.out rerun_struct11_coord6.xvg rerun_struct3_coord7.xvg rerun_struct6_coord8.xvg rerun_struct9_coord9.xvg rerun_struct11_coord7.edr rerun_struct3_coord8.edr rerun_struct6_coord9.edr unk_matrix.pkl rerun_struct11_coord7.log rerun_struct3_coord8.log rerun_struct6_coord9.log rerun_struct11_coord7.out rerun_struct3_coord8.out rerun_struct6_coord9.out rerun_struct11_coord7.xvg rerun_struct3_coord8.xvg rerun_struct6_coord9.xvg

(fep) cadd@huiyu:~/program/fep/workdir/fep$ bash runall.sh

environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
environment: rerun/rerun_struct_coord.out: No such file or directory
tar: /water/md/lambda/analysis.log: Cannot stat: No such file or directory
tar: /protein/md/lambda/analysis.log: Cannot stat: No such file or directory
tar: /water/md/analysis/.xvg: Cannot stat: No such file or directory
tar: /protein/md/analysis/.xvg: Cannot stat: No such file or directory
tar: Exiting with failure status due to previous errors

Originally posted by @Lumen97 in #56 (comment)

(PyAutoFEP) tbaba@YW113:/work113/test/PyAutoFEP/docs/tutorial01_fin$ python ../../analyze_results.py --input tutorial.tgz --units kcal --output_uncompress_directory tutorial --center_molecule=FXR_12

Bad key "text.kerning_factor" on line 4 in
/home1/test/ubuntu20/anaconda3/envs/PyAutoFEP/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test_patch.mplstyle.
You probably need to get an updated matplotlibrc file from
http://github.com/matplotlib/matplotlib/blob/master/matplotlibrc.template
or from the matplotlib source distribution
All available analysis will be run
Traceback (most recent call last):
File "../../analyze_results.py", line 1478, in
verbosity=arguments.verbose)
File "../../analyze_results.py", line 1139, in collect_from_xvg
with open(os.path.join(xvg_directory, unk_file), 'rb') as fh:
NotADirectoryError: [Errno 20] Not a directory: 'tutorial/FXR_12-FXR_74/runall_FXR_12-FXR_74_water.sh/md/rerun/unk_matrix.pkl'

After all the calculations were completed, I ran the above command and got an error.
The cause is the path of the file.　It should be tutorial/FXR_12-FXR_74/water/md/rerun/unk_matrix.pkl, but it is tutorial/FXR_12-FXR_74/runall_FXR_12-FXR_74_water.sh/md/rerun/unk_ matrix.pkl. I was not sure where to modify to change the path name, so I would appreciate your advice.

Hi, it is mentioned in the manual ‘It automates the building of a simple receptor-ligand system and will fail for complex cases (eg: membrane systems, systems with vacuum regions, solvent mixtures, and systems with nucleic acids on cofactors).’ So I want to know whether PyAutoFEP can be applied to the calculation of membrane system. Looking forward to your reply, thanks!

Hi there, thank you for maintaining this awesome repo! I'm wondering whether it is possible to calculate absolute binding free energy with PyAutoFEP, either with the current code base or some modifications on it.

When using matplotlib 2.2.3 py37hb69df0a_0, analyze_results.py fails with

All available analysis will be run
================== Pairwise ΔΔG ==================
         Perturbation           protein    water  
[...]/anaconda3/envs/pyautofep_temp/lib/python3.7/site-packages/matplotlib/tight_layout.py:228: UserWarning: tight_layout : falling back to Agg renderer
  warnings.warn("tight_layout : falling back to Agg renderer")
Traceback (most recent call last):
  File "[...]/PyAutoFEP/analyze_results.py", line 1645, in <module>
    ddg_data_result = [analyze_perturbation(**each_kwargs) for each_kwargs in arguments_data_holder.values()]
  File "[...]/PyAutoFEP/analyze_results.py", line 1645, in <listcomp>
    ddg_data_result = [analyze_perturbation(**each_kwargs) for each_kwargs in arguments_data_holder.values()]
  File "[...]/PyAutoFEP/analyze_results.py", line 935, in analyze_perturbation
    verbosity=verbosity)
  File "[...]/PyAutoFEP/analyze_results.py", line 788, in plot_stacked_bars
    ticks=numpy.arange(0, n_rep) + 0.5)
  File "[...]/anaconda3/envs/pyautofep_temp/lib/python3.7/site-packages/matplotlib/pyplot.py", line 2320, in colorbar
    ret = gcf().colorbar(mappable, cax = cax, ax=ax, **kw)
  File "[...]/anaconda3/envs/pyautofep_temp/lib/python3.7/site-packages/matplotlib/figure.py", line 2098, in colorbar
    cb = cbar.colorbar_factory(cax, mappable, **cb_kw)
  File "[...]/anaconda3/envs/pyautofep_temp/lib/python3.7/site-packages/matplotlib/colorbar.py", line 1399, in colorbar_factory
    cb = Colorbar(cax, mappable, **kwargs)
  File "[...]/anaconda3/envs/pyautofep_temp/lib/python3.7/site-packages/matplotlib/colorbar.py", line 920, in __init__
    mappable.autoscale_None()
  File "[...]/anaconda3/envs/pyautofep_temp/lib/python3.7/site-packages/matplotlib/cm.py", line 362, in autoscale_None
    raise TypeError('You must first set_array for mappable')
TypeError: You must first set_array for mappable

It works with matplotlib 3.4.2 py37h06a4308_0.

Dear Prof.

We finished the FEP calculations of FXR (12 to 74, 76, 84, 85, 88) without REST2, but got different results from you(tutorial01/tutorial/ddg_to_center.csv). We also noticed that the calculated results of 5 ligands are different from the experimental values in the article.
Our calculation results are shown as follows：

  | Our ddg | Tutorial | Exp. Dg
FXR_12 | 0 | 0 | -9.93271
FXR_74 | 0.070849 | 0.109407 | -8.48765
FXR_76 | 0.185977 | -0.23893 | -6.00989
FXR_84 | 0.009766 | -0.16105 | -7.33357
FXR_85 | 0.079347 | 0.511437 | -8.95884
FXR_88 | 0.272636 | 0.759893 | -8.60273

For FEP calculations, all parameters and procedures we use are based on those given in tutorial [OPLS/AA-m, Ligpargen, gromacs, plumed, without REST2/hrex]
Our .ini configuration is shown below

# This is the configuration section for prepare_dual_topology.py (which is the only section in this file, but the [ section ] is mandatory)
[prepare_dual_topology]

# Read ligands and topologies from this folder
input_ligands = lig_data

# Read the macromolecule structure from this file
structure = receptor_data/5q17_processed.pdb

# This is the force field directory, it will be copied to each perturbation dir and used to prepare the MD systems
extradirs = oplsaam.ff

# Options controlling the core-constrained superimposition
# First select the use of it instead of reading all ligand poses
pose_loader = superimpose

#lambda set up
lambda_input = lambdas12

# Use this pose as the reference for the superimposition
poses_reference_pose_superimpose = receptor_data/9mv.pdb

# Name of the output. This will be a self-extracting bash file
perturbations_dir = tutorial-t1

# Sets the path to GROMACS executable in the run node, uncomment and modify if needed.
# gmx_bin_run = /usr/local/bin/gromacs

output_hidden_temp_dir=False

# Options controlling the output, see manual for more info. Uncomment as needed
# FEP legs are to be submitted to a slurm scheduler
#output_scripttype = slurm
# Run these commands at the beginnig of the jog (useful to load modules, importing libs)
# output_runbefore = module load python3; module load cuda
# Fine-tune job resources
#output_resources = all_cpus:48; all_gpus:4; all_time: 24
# Use a python file instead of a binary during the collect step
# output_collecttype = python

So we would like to ask you if you can give some guidance on how to calculate the simulation accurately. And qe would be deeply grateful if you could give more detailed configuration files of this case.

In case the MCS matches more than a set of atoms in one (or in both) the ligands, it is possible that the atom matching between MCS and the ligand will be different than the on that was matched during the calculation of the MCS. This will lead to a incorrect dual-topology. I believe that the only solution to this problem is to check whether each of the matches would yield a correct dual-topology. In case there are multiple matches for both ligands, a combinatorial problem arises.

After the last line above, it would be better check if xvg_directory is in valid status like:

            xvg_directory = os.path.join(data_dir, each_entry, system, 'md', 'rerun')
            if not os.path.isdir(xvg_directory) or len(glob(os.path.join(xvg_directory,'*.xvg'))) == 0:
                continue

In case the directory is not a directory or xvg files are not found, it can be skipped. The first case occurs if you run analyze_results.py without packing the output directory,. The runall*.sh files in the leg directory give errors without the modification above.

The software "parallel" should be mentioned in the installation part.

Further, I use a single node workstation. I need to modify some commands. I'm not sure if it's because of my software environment.
For runall_protein.sh, line 221, "gmx mdrun -v" need to be changed to "gmx_mpi mdrun -v"
For runall_water.sh, line 197, "gmx mdrun -v" need to be changed to "gmx_mpi mdrun -v"

Also, it would be great to couple that with the automatic lambda scheme selection.

If I choose the kBT unit, the beta is 1/temperature, but the unit of the mbar estimator_obj delta_f_ is already kT.
why output ddg = delta_f_/(1/temperature) = delta_f_*temperature

Hi all,

while running this

<sbatch --parsable FXR_12-FXR_74/runall_FXR_12-FXR_74_protein.sh>

the following error showed up, anyone come-across this error before or any suggestion?

Thank you!

GROMACS: gmx mdrun, version 2021.5
Executable: /home/easybuild/.local/easybuild/software/GROMACS/2021.5-foss-2021b-CUDA-11.4.1/bin/gmx_mpi
Data prefix: /home/easybuild/.local/easybuild/software/GROMACS/2021.5-foss-2021b-CUDA-11.4.1
Working dir: /home/gdegaga_championsoncology_com/PyAutoFEP/Slurmtest/tutorial/slurm00/FXR_12-FXR_74/protein/md
Command line:
gmx_mpi mdrun -v --deffnm lambda -cpi lambda.cpt -multidir lambda0 lambda1 lambda2 lambda3 lambda4 lambda5 lambda6 lambda7 lambda8 lambda9 lambda10 lambda11

Program: gmx mdrun, version 2021.5

Program: gmx mdrun, version 2021.5
Source file: src/gromacs/mdrunutility/handlerestart.cpp (line 575)
Function: void gmx::{anonymous}::StartingBehaviorHandler::ensureMultiSimBehaviorsMatch(const gmx_multisim_t*)
MPI rank: 1 (out of 12)

Inconsistency in user input:

Multi-simulations must all start in the same way, either a new
simulation, a restart with appending, or a restart without appending.
However, the contents of the multi-simulation directories you specified
were inconsistent with each other. Either remove the checkpoint file
from each directory, or ensure each directory has a checkpoint file from

I'm currently creating input for another system.
I have a few questions about Ligand and I would like to share them with you.

I have found that some input files can be created with the same structure, but some cannot, depending on the order of the atoms. Is there anything I should pay attention to when preparing the files?
　　Success LIG13.mol -> LIG0.mol
　　Failure LIG13_new.mol→LIG0_new.mol

https://github.com/baba-hashimoto/PyAutoFEP/tree/master/files
[ERROR] Less than 3 common atoms were found between molecule LIG13 (SMILES=<rdkit.Chem.rdchem.Mol object at 0x14562bb603a0>) and LIG0 (SMILES=(17, 17))

More than half of the constituent atoms belong to dual_topology, have you ever calculated such a system? In this case, I was able to prepare the input file, but I got an error in the energy minimization calculation.

[ exclusions ] ; Added by dual_topology
11 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 ; Added by dual topology
12 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 ; Added by dual topology

I'm interested in reading the mapping and using it in my analysis. I wonder if the PyAutoFEP could generate such a mapping?
For example

Molecule 1 : 0 1 2 3
               | | |
Molecule 2:    0 1 2 3

Would generate such a mapping file

mapping = {('Molecule_1', 'Molecule_2'): [(0, None),
                                          (1, 0),
                                          (2, 1),
                                          (3, 2),
                                          (None, 3),]}

Now that the without rest2 calculation works, I am trying the with rest2 calculation.
https://www.plumed.org/doc-v2.7/user-doc/html/hrex.html　
I understand that Tools/rest2.sh is running to perform the task described in the link above.

Currently, I am getting the following error.

[ERROR] Failed to read scaling selection file.
[ERROR] Could not read file LIG: O00_const_idx1 (and read_file_to_buffer wall called with die_on_error=False). Error was: [Errno 2] No such file or directory: 'LIG: O00_const_idx1'
=================== STACK INFO ===================
File "../../prepare_dual_topology.py", line 2948, in
solute_scaling_atoms_dict = process_scaling_input(arguments.solute_scaling_selection)
File "../../prepare_dual_topology.py", line 2160, in process_scaling_input
'', verbosity=verbosity)
File "/work113//PyAutoFEP/os_util.py", line 99, in read_file_to_buffer
local_print(error_message, msg_verbosity=verbosity_level.error, current_verbosity=verbosity)
File "/work113//PyAutoFEP/os_util.py", line 292, in local_print
formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================

As far as the program is concerned, it fails to write solute_scaling_slection. If you could give me an example of how to fill it out, that would be helpful.

Hi! I read your journal named PyAutoFEP: An Automated Free Energy Perturbation Workflow for GROMACS Integrating Enhanced Sampling Methods, which is such a great work.

I try to do some FEP calculations in your method while i met some problems.
I used the smis.smi(see below) as input of the generate_perturbation_map.py with command line.

python generate_perturbation_map.py --input smis.smi

the smis.smi is:

CN1C(=NC2(CC(C)(C)Oc3ccc(cc23)c4cccc(Cl)c4)C1=O)N aname
CN1C(=NC2(CC(C)(C)Oc3cnc(cc23)c4cccc(Cl)c4)C1=O)N bname
CN1C(=NC2(CC(C)(C)Oc3ncc(cc23)c4cccc(Cl)c4)C1=O)N cname
CN1C(=NC2(CC(C)(C)Oc3ccc(cc23)c4cccc(I)c4)C1=O)N dname
CN1C(=NC2(CC(C)(C)Oc3ccc(cc23)c4cccc(Br)c4)C1=O)N ename

However, i met a TypeError of :

Traceback (most recent call last):
  File "generate_perturbation_map.py", line 497, in <module>
    custom_mcs=custom_user_data, savestate=progress_data, verbosity=arguments.verbose)
  File "generate_perturbation_map.py", line 138, in fill_thermograph
    for (edge_i, edge_j) in thermograph.edges:
TypeError: 'method' object is not iterable

Did I use that in a wrong way?
Thank you for your help.

I'm interested in running the tutorial (https://github.com/luancarvalhomartins/PyAutoFEP/tree/master/docs/tutorial01).
I have executed the ~/GitHub/PyAutoFEP/generate_perturbation_map.py --map_type=star --map_bias=FXR_12 --input lig_data/*.mol, which runs fine.
Then I run ~/GitHub/PyAutoFEP/prepare_dual_topology.py --config_file=step2.ini and I got

=============== mdp and run steps=================
         Complex                    Water          
        min01.mdp                 min01.mdp        
        min02.mdp                 min02.mdp        
        min03.mdp                  nve.mdp         
         nve.mdp                   nvt.mdp         
         nvt.mdp                   npt.mdp         
         npt.mdp                   md.mdp          
         md.mdp                       -            
==================================================
================= Input ligands ==================
   Name         Molecule            Topology      
  FXR_12   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Cl)c1[H] lig_data/FXR_12.itp
  FXR_74   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(Br)c1[H] lig_data/FXR_74.itp
  FXR_76   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c([H])c1[H] lig_data/FXR_76.itp
  FXR_84   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(F)c1[H] lig_data/FXR_84.itp
  FXR_85   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C([H])([H])[H])c1[H] lig_data/FXR_85.itp
  FXR_88   [H]c1c([H])c([H])c(S(=O)(=O)N2C([H])([H])C([H])([H])C3(C(=O)N(C([H])([H])c4c([H])c([H])c(C(=O)[O-])c([H])c4[H])c4c([H])c([H])c(Br)c([H])c43)C([H])([H])C2([H])[H])c(C(F)(F)F)c1[H] lig_data/FXR_88.itp
==================================================
=================== Poses read ===================
Name            File                
reference       {'molecule': 'receptor_data/9mv.pdb', 'format': '.pdb'}
=================== Poses read ===================
Name            File                
FXR_12          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eac10>, 'topology': ['lig_data/FXR_12.itp']}
FXR_74          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207ead00>, 'topology': ['lig_data/FXR_74.itp']}
FXR_76          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eabc0>, 'topology': ['lig_data/FXR_76.itp']}
FXR_84          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eacb0>, 'topology': ['lig_data/FXR_84.itp']}
FXR_85          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207ead50>, 'topology': ['lig_data/FXR_85.itp']}
FXR_88          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eac60>, 'topology': ['lig_data/FXR_88.itp']}
=============== Superimposed poses ===============
Name            File                      Note           
FXR_12          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eac10>, 'topology': ['lig_data/FXR_12.itp']}
FXR_74          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207ead00>, 'topology': ['lig_data/FXR_74.itp']}
FXR_76          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eabc0>, 'topology': ['lig_data/FXR_76.itp']}
FXR_84          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eacb0>, 'topology': ['lig_data/FXR_84.itp']}
FXR_85          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207ead50>, 'topology': ['lig_data/FXR_85.itp']}
FXR_88          {'molecule': <rdkit.Chem.rdchem.Mol object at 0x7fe3207eac60>, 'topology': ['lig_data/FXR_88.itp']}
==================================================
================== Align poses ===================
Molecule FXR_12 aligned
Molecule FXR_74 aligned
Molecule FXR_76 aligned
Molecule FXR_84 aligned
Molecule FXR_85 aligned
Molecule FXR_88 aligned
================= Perturbations ==================
         State A                  State B         
         FXR_12                   FXR_74          
         FXR_12                   FXR_76          
         FXR_12                   FXR_84          
         FXR_12                   FXR_85          
         FXR_12                   FXR_88          
==================================================
================ Working on pairs ================
Perturbation Pose Coordinates
============ Working on FXR_12-FXR_74 ============
FXR_12 -> FXR_74
================ Building system =================
[ERROR] Failed to run gmx grompp. Error code 1.
[ERROR] Command line was: ['gmx', 'grompp', '-f', '/var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/genion_step5_123555_18102021.mdp', '-c', '/var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/systemsolvated_step4_123555_18102021.gro', '-p', '/var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/systemsolv_step4_123555_18102021.top', '-o', '/var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/genion_step6_123555_18102021.tpr', '-maxwarn', '1', '-po', '/var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/mdout_step6_123555_18102021.mdp']
[ERROR] 
[ERROR] stdout:
[ERROR] Setting the LD random seed to 2009067503
[ERROR] 
[ERROR] Generated 449826 of the 449826 non-bonded parameter combinations
[ERROR] 
[ERROR] Generated 449826 of the 449826 1-4 parameter combinations
[ERROR] 
[ERROR] Excluding 3 bonded neighbours molecule type 'Protein'
[ERROR] 
[ERROR] Excluding 3 bonded neighbours molecule type 'Protein2'
[ERROR] 
[ERROR] Excluding 3 bonded neighbours molecule type 'LIG'
[ERROR] 
[ERROR] 
[ERROR] stderr:
[ERROR]                       :-) GROMACS - gmx grompp, 2021.3 (-:
[ERROR] 
[ERROR]                             GROMACS is written by:
[ERROR]      Andrey Alekseenko              Emile Apol              Rossen Apostolov     
[ERROR]          Paul Bauer           Herman J.C. Berendsen           Par Bjelkmar       
[ERROR]        Christian Blau           Viacheslav Bolnykh             Kevin Boyd        
[ERROR]      Aldert van Buuren           Rudi van Drunen             Anton Feenstra      
[ERROR]     Gilles Gouaillardet             Alan Gray               Gerrit Groenhof      
[ERROR]        Anca Hamuraru            Vincent Hindriksen          M. Eric Irrgang      
[ERROR]       Aleksei Iupinov           Christoph Junghans             Joe Jordan        
[ERROR]     Dimitrios Karkoulis            Peter Kasson                Jiri Kraus        
[ERROR]       Carsten Kutzner              Per Larsson              Justin A. Lemkul     
[ERROR]        Viveca Lindahl            Magnus Lundborg             Erik Marklund       
[ERROR]         Pascal Merz             Pieter Meulenhoff            Teemu Murtola       
[ERROR]         Szilard Pall               Sander Pronk              Roland Schulz       
[ERROR]        Michael Shirts            Alexey Shvetsov             Alfons Sijbers      
[ERROR]        Peter Tieleman              Jon Vincent              Teemu Virolainen     
[ERROR]      Christian Wennberg            Maarten Wolf              Artem Zhmurov       
[ERROR]                            and the project leaders:
[ERROR]         Mark Abraham, Berk Hess, Erik Lindahl, and David van der Spoel
[ERROR] 
[ERROR] Copyright (c) 1991-2000, University of Groningen, The Netherlands.
[ERROR] Copyright (c) 2001-2019, The GROMACS development team at
[ERROR] Uppsala University, Stockholm University and
[ERROR] the Royal Institute of Technology, Sweden.
[ERROR] check out http://www.gromacs.org for more information.
[ERROR] 
[ERROR] GROMACS is free software; you can redistribute it and/or modify it
[ERROR] under the terms of the GNU Lesser General Public License
[ERROR] as published by the Free Software Foundation; either version 2.1
[ERROR] of the License, or (at your option) any later version.
[ERROR] 
[ERROR] GROMACS:      gmx grompp, version 2021.3
[ERROR] Executable:   /usr/local/gromacs/2021.3/bin/gmx
[ERROR] Data prefix:  /usr/local/gromacs/2021.3
[ERROR] Working dir:  /Volumes/GoogleDrive/My Drive/nutshell/workdir
[ERROR] Command line:
[ERROR]   gmx grompp -f /var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/genion_step5_123555_18102021.mdp -c /var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/systemsolvated_step4_123555_18102021.gro -p /var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/systemsolv_step4_123555_18102021.top -o /var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/genion_step6_123555_18102021.tpr -maxwarn 1 -po /var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/mdout_step6_123555_18102021.mdp
[ERROR] 
[ERROR] 
[ERROR] NOTE 1 [file /var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/genion_step5_123555_18102021.mdp]:
[ERROR]   For a correct single-point energy evaluation with nsteps = 0, use
[ERROR]   continuation = yes to avoid constraining the input coordinates.
[ERROR] 
[ERROR] Generating 1-4 interactions: fudge = 0.5
[ERROR] 
[ERROR] NOTE 2 [file systemsolv_step4_123555_18102021.top, line 56]:
[ERROR]   System has non-zero total charge: -10.000000
[ERROR]   Total charge should normally be an integer. See
[ERROR]   http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
[ERROR]   for discussion on how close it should be to an integer.
[ERROR]   
[ERROR] 
[ERROR] 
[ERROR] 
[ERROR] -------------------------------------------------------
[ERROR] Program:     gmx grompp, version 2021.3
[ERROR] Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 681)
[ERROR] 
[ERROR] Fatal error:
[ERROR] number of coordinates in coordinate file
[ERROR] (/var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/systemsolvated_step4_123555_18102021.gro,
[ERROR] 44742)
[ERROR]              does not match topology
[ERROR] (/var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/systemsolv_step4_123555_18102021.top,
[ERROR] 3930)
[ERROR] 
[ERROR] For more information and tips for troubleshooting, please check the GROMACS
[ERROR] website at http://www.gromacs.org/Documentation/Errors
[ERROR] -------------------------------------------------------
[ERROR] 
[ERROR] 
[ERROR] ==================================================
[ERROR] This is likely caused by a failing to edit the intermediate topology file /var/folders/31/qmmnwsvd13j74h338x3d57w8000d3z/T/tmpmye3fjzm/tutorial/FXR_12-FXR_74/protein/build_system_123555_18102021/systemsolv_step4_123555_18102021.top. Rerunning with output_hidden_temp_dir=False may solve this issue.
=================== STACK INFO ===================
  File "/Users/grte2001/GitHub/PyAutoFEP/prepare_dual_topology.py", line 3290, in <module>
    solvate_data=solvate_data, verbosity=arguments.verbose)
  File "/Users/grte2001/GitHub/PyAutoFEP/prepare_dual_topology.py", line 480, in prepare_complex_system
    msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
  File "/Users/grte2001/GitHub/PyAutoFEP/os_util.py", line 292, in local_print
    formatted_string += '\n{:=^50}\n{}{:=^50}'.format(' STACK INFO ', ''.join(traceback.format_stack()),
=================== STACK INFO ===================

I wonder if I can get some help with regard to this, please? Thank you.

Hi, I encoutered Could not find topology data for ligand xxx.pdb error at the prepare_dual_topology.py --config_file=step2.ini stage. Detailed massages are in stdout.log within the attachment FEP.zip. I think I have well prepared all the topology files, which are in lig_data dir

[FEP.zip](https://github.com/luancarvalhomartins/PyAutoFEP/files/9394041/FEP.zip)

During the rerun step, gmx mdrun -rerun is used to obtain potential energies of each of the sampled conformations of the system using every hamiltonian. During this process, if the potential energy is too high, an error like the one below (note: energies will vary) will be raised and gmx will halt execution.

Fatal error:
Step 0: The total potential energy is 8.11303e+09, which is extremely high.
The LJ and electrostatic contributions to the energy are 8.11306e+09 and
-31622, respectively. A very high potential energy can be caused by
overlapping interactions in bonded interactions or very large coordinate
values. Usually this is caused by a badly- or non-equilibrated initial
configuration, incorrect interactions or parameters in the topology.

This error, is, therefore, a showstopper. gmx mdrun -rerun only checks for this high potential in the first frame, so a workaround could be removing the first frame of the trajectories and retrying, which feels a bit hacky. So far, I found no other solution, though. In case anyone have any suggestion about this, feel free to comment.

Note: in case you are affected by this issue and a fix is yet to be committed, use the following steps to automatically prepare trajs without the first frame and rerun gmx mdrun -rerun, if needed.

> parallel --jobs=1 'printf "0\n" | gmx trjconv -f {}/lambda.trr -s {}/lambda_nocont.tpr -o {}/lambda_skip0.xtc -b 0.002 -ndec 6' ::: lambda*
> parallel --jobs=1 'gmx dump -e rerun/rerun_struct{2}_coord{1}.edr &> /dev/null || gmx mdrun -rerun lambda{1}/lambda_skip0.xtc -s lambda{2}/lambda_nocont.tpr -e rerun/rerun_struct{2}_coord{1}.edr -g rerun/rerun_struct{2}_coord{1}.log' ::: $(seq 0 11) ::: $(seq 0 11)

This error showed up suddenly, any reason and solution to this would be greatly appreciated.

PyAutoFEP]$ python prepare_dual_topology.py --config_file step2.ini
Traceback (most recent call last):
File "/home/gdegaga_championsoncology_com/PyAutoFEP/prepare_dual_topology.py", line 2952, in
pybel.ob.obErrorLog.SetOutputLevel(pybel.ob.obError)
AttributeError: module 'pybel' has no attribute 'ob'

Thanks!

Hi,
First of all, thank you for sharing your precious source code.

I am having trouble in topology generation regarding solvation.

Please let me know if the issues can be resolved without changing the source code.

Checklist for closing this (Added by @luancarvalhomartins so GitHub can track completion, hope you don't mind me editing your issue):

• Allowing for setting ionic strength of the water leg (independent of complex leg)
• Passing arguments.buildsys_water to guess_water_box (1074ad4)
• Allow for arbitrary water molecule names (Added 2022-02-09) (189fc2f)
• Run gmx genion for solvated systems
• Rewrite parts of prepare_dual_topolopy.py to avoid assuming the molecule name of the protein (how exactly?)

With the new release, the issue (choderalab/pymbar#419) is resolved by the more robust MBAR solver, which could be called as

from alchemlyb.estimators import AutoMBAR as MBAR

Dear Developer

Q1: I am using the analyze_results.py to evaluate the delta detla G (DDG). While learning and debugging the code, I notice there may be an issue while collecting DDG after calling " estimator_obj.fit(subsampled_u_nk_after_first)" as below

this_ddg = -1 * estimator_obj.delta_f_[0].iloc[-1] / beta

In my case, the number of total states ( including the start, end and intermedium states ) is 12 ( index: 0 - 11 ). Here is the delta_f_ in my case.

When I check the return value of estimator_obj.delta_f_[0].iloc[-1], it gives

The value corresponds to the last row of the second column (in red). I think this is DDG between state 0 and state 7. It seems the right one should be DDG = -4.667543 between state 0 and state 11 (the fifth row of the second column, in blue). The problem is the states is not in the ascending order.

Could you kindly help me to check if there is a problem?

Q2: I also like to report the other issue. I found a logical issue in the following code.

It seems the code will always run the two blocks in blue.

Thanks,
ZR

According to some energy minimization and MD mdp files, positional restrains are desirable. However, the topology files do not have [ position_restraints ] in both complex and ligand systems. Also, the definition has a different name. In the complex topology generated, the definition name is POSRES not POSRES_PROTEIN. Input files such as min01.mdp requires POSRES_PROTEIN as the definition name.
I manually add [ position_restraints ] after calling prepare_dual_topology, but it looks like prepare_dual_topology should handel it.

Thank very much for sharing your great work. I've got a question when reading the codes and hope you can explain the logic behind.

When embed it to the reference pose, constrained_embed_dualmol() function was used to generate an embedding of a dual-topology pseudomolecule to target using MCS. You commented 'Prepare a temporary copy of molecules A and B, but remove original conformers to make sure we don't get one of them in the rms_dict' .

My question is: Why remove all the conformers of molecule A and B and replace them with newly generated conformers?

In my case, the input target is often molecule A itself and has already be aligned in align_ligands function during superimpose step. Therefore, keep the superimposed poses_mol_data would be a better option for me, especially num_conformers=50 used here is less than num_conformers=200 used in align_ligands step, which result in worse alignment with higher volume function values sometimes.

I found a mistake in output_files_data.ini.
I think line 199 should be changed as follows.
( parallel --jobs=$n_jobs 'GMXBIN distance -select '''com of group LIG plus com of group Protein''' -f {.} _pbcfit.xtc -s {} -n ... /INDEX -oall analysis/{//}_protligdist.xvg >> {//}/analysis.log 2>&1' ::: ${trajectory_files} ) || { echo "Failed to calculate distance between ligand and protein on line $LINENO"; }
Could you check the script ?

PyAutoFEP dies with an error when output_scripttype = pbs is used.

================ Working on pairs ================
Perturbation Pose Coordinates
Traceback (most recent call last):
  File "prepare_dual_topology.py", line 3189, in <module>
    verbosity=arguments.verbose)
  File "prepare_dual_topology.py", line 675, in prepare_output_scripts_data
    temp_str = template_section['header']
  File "/home/luancarvalho/anaconda3/envs/PyAutoFEP_testrepo/lib/python3.7/configparser.py", line 1251, in __getitem__
    raise KeyError(key)
KeyError: 'header'