find_mcs function could include different atom types with the same enviroment into the resulted MCS

[YunpengCai]

Hi, I faced some failure during MD simulation with error that one or more water molecules cannot be settled. It could be the system I prepared is somehow incorrect.

After looking into the dual molecule prepared in the system, I think it might improve the system if we get the right MCS from a CADD scientist's view. For instance, ligand_A O=C(NC1=CC2=C(CC(CC2)C2CCC2)C=C1)C1CC1 and ligand_B FC(F)(F)C(=O)NC1=CC2=C(CC(NC2)C2CCC2)C=C1 will not find the largest MCS due to the difference of C and N atoms on the ring.

I am not sure whether it is the reason that fails my simulation and would ask about your opinion about how the MCS will affect the results. Thanks a lot.

[luancarvalhomartins]

First, is your MCS O=CNc1ccccc1? In case it isn't, let me know. Maybe just fixing this could solve your problem.

It is quite likely that the size of the common region is too small and that is leading to instabilities on the simulation. However, it is also possible that this is not the reason and the problem is something else (eg, the starting pose is not stable enough, the system requires more minimization and/or equilibration). It's not always simple to know.

Currently, PyAutoFEP cannot (a) perturb bonds and (b) concomitantly perturb charges and VdW parameters of an atom; so that it is not possible to perturb a C to a N in ring right now. That's why the MCS code won't allow for such core (it would give spurious, unphysical results). I am finishing a code for doing charge-changing perturbations and that code will allow for such perturbations to be done (as a side effect of the implementation). Should you be able to wait for a few weeks, it will make into the tree.

[luancarvalhomartins]

Allow for perturbing atoms in rings

• [x] Code for concomitant perturbation of charges and VdW (done in charged perturbations branch)
• [ ] Code for perturbing bonds (done in charged perturbations branch)
• [ ] Update MCS code
• [ ] Validate

[luancarvalhomartins]

(It seems you deleted your message, so I don't know if this is still relevant. Nevertheless, please find below my answers.)

Doing position restraints can improve the stability of the system, however some caution is required when applying it. It will prevent FEP from exploring conformational changes of the ligand. Applying restraints on the A or B part of the ligand can have unpredictable effects, as PyAutoFEP won't correct for any change in the number of degrees of freedom introduced by the restraints.

By the design of GROMACS topology/structure implementation, solute_scaling_selection must match the atoms in the topology. You can match both ligand and protein atoms (but not the solvent ones). Some examples:

# Apply scaling to all ligand atoms
solute_scaling_selection = LIG: *

# Apply scaling to some ligand atoms only using atom indexes (pass a list, must match the number in the topology, starts form 1)
solute_scaling_selection = LIG: 1,2,3,4,5,10,11,12,20,25

#  Apply scaling to some ligand atoms only using atom types (must match the names in the topology, regex allowed)
solute_scaling_selection = LIG: C0K_const_idx21, C0M_const_idx22, C0N_const_idx23

#  Apply scaling to some ligand atoms using atom types and regex (must match the names in the topology)
solute_scaling_selection = LIG: *_top*

#  Apply scaling to some protein atoms (avoid using atom types as they would atoms in any res)
solute_scaling_selection = Protein: 2943, 2944, 2945, 2946, 2947, 2948, 2949, 2950, 2951, 2952, 2953, 2954, 2955, 2956, 2957, 2958, 2959

The documentation on solute scaling should be improved. Hope I can address it atsome point.