Force fields

BindFlow provides several out-of-the-box force fields. This section explains how to access them and integrate custom force fields within BindFlow. We mainly focus on the input possibilities for BindFlow’s runner functions.

BindFlow offers a variety of force field options, but as Uncle Ben says, “With great power comes great responsibility.” Users must choose the appropriate combination of force fields. By default, BindFlow offers a suitable combination.

Structure inputs

Six keywords control the type of force field used for each specific component in the system:

1. protein: Definition of the host

2. membrane: Definition of the membrane

3. ligands: A list of ligand’s definitions

4. cofactor: Definition of the cofactor

5. water_model: Type of GROMACS’ water model to use

6. custom_ff_path: Path to the custom force field if needed.

Partial and Full Definitions

For a straightforward setup, you can provide the path to the corresponding file(s), which we will call the partial definition. However, you also have the option to fine-tune the definition of force fields for each component, referred to as the full definition.

The quality of the initial structure is critical for accurate results, whether using FEP or MM(PB/GB)SA Behera et al, 2025. Proper definitions of tautomeric, isomeric, and protonation states for both proteins and ligands are essential as well the conformation state of the protein-ligand complex. BindFlow does not aim to solve this issue directly, as specialized tools excel in this domain. However, BindFlow offers basic functionality to “fix” proteins (e.g., resolving missing atoms or correcting atom naming) using pdbfixer from OpenMM and pdb2gmx from GROMACS. For complex receptors, users are advised to provide fully defined structures (e.g., GRO and TOP files) or preprocessed, compatible PDB files.

In the following examples, we will use the runner bindflow.runners.calculate().

Partial definition

protein

The protein will be processed with amber99sb-ildn force field after been fixed (if fix_protein = True) with PDBFixer.

Hint

It is advised to spend some time on the processing of the protein beforehand (better a minute than repeat the whole campaign):

1. Missing atoms

2. Missing loops

3. Terminal capping

4. Protonation state

All the above steps are highly system-dependent, and while PDBFixer can handle some minor issues, it is far from perfect. In addition, our use of PDBFixer is very simple.

calculate(
    ...
    protein="path/to/protein.{pdb;gro}",
    ...   
)

membrane

The membrane will be processed with SLipids_2020.

Getting the membrane.pdb file

For a membrane systems you first need to embed the protein into the membrane. This can easily done with CHARMM-GUI. This is a non-exhaustive list of steps for this process:

Processing on CHARMM-GUI:

• ACE and CT3 terminus

• pH=7

• Run PPM 2.0

• Hexagonal box

• Water thickness: 15

• Length of X and Y: 90 (initial guess)

• Only POPC for a simple membrane or any other lipid composition you would like to use

• Do not include ions

• force field:

  • AMBER:

    • Protein: ff14sb

    • Lipid: SLipids

    • Water: TIP3P

• Input Generation Options: GROMACS

Open gromac/step5_input.gro. This file has the crystal information which is important during the solvation step in BindFLow. Convert to PDB with gmx editconf -f step5_input.gro -o temporal.pdb . The last PDB will also have the crystal info. Split the PDB in POPC and protein.

In PyMOL:

• select popc, resn POPC

• select prot, (polymer.protein or resn ACE or resn NME)

👀 Depending on how you are processing the protein, you might need to change the entry ATOM to HETATM for the CAP groups ACE and NME manually in the PDB file.

calculate(
    ...
    membrane="path/to/membrane.pdb",
    ...   
)

ligands

calculate(
    ...
    ligands=[
        "path/to/ligand1.{mol;sdf}",
        "path/to/ligand2.{mol;sdf}",
        "path/to/ligand3.{mol;sdf}",
        ...
    ],
    ...   
)

cofactor

calculate(
    ...
    cofactor="path/to/cofactor.{mol;sdf}",
    ...   
)

Full definition

protein

by code

on GROMACS distribution

You can access all the GROMACS force fields by their code, they will be pass to pdb2gmx through the flag -ff after been fixed (if fix_protein = True) with PDBFixer.

Hint

It is advised to spend some time on the processing of the protein beforehand (better a minute than repeat the whole campaign):

1. Missing atoms

2. Missing loops

3. Terminal capping

4. Protonation state

All the above steps are highly system-dependent, and while PDBFixer can handle some minor issues, it is far from perfect. In addition, our use of PDBFixer is very simple

calculate(
    ...
    protein={
        "conf": "path/to/protein.{pdb;gro}",
        "ff":{
            "code": <GMX_force_field_code>
        }
    },
    ...   
)

external

To add even more flexibility, you can use any external force field ported to GROMACS, in this case you just need to copy your force_field.ff (e.g. charmm36-jul2022.ff) to your desired directory and pass the path of this directory to custom_ff_path parameter. If you have more force fields, you can copy all of them in the same directory. BindFlow will internally set the following environmental variable at run time.

os.environ["GMXLIB"] = os.path.abspath(custom_ff_path)

Warning

See how the force field directory ends in .ff; e.g. charmm36-jul2022.ff. This is needed.

The force field code (e.g. for charmm36-jul2022.ff, the code is charmm36-jul2022) will be pass to pdb2gmx through the flag -ff after been fixed (if fix_protein = True) with PDBFixer.

Hint

It is advised to spend some time on the processing of the protein beforehand (better a minute than repeat the whole campaign):

1. Missing atoms

2. Missing loops

3. Terminal capping

4. Protonation state

All the above steps are highly system-dependent, and while PDBFixer can handle some minor issues, it is far from perfect. In addition, our use of PDBFixer is very simple

calculate(
    ...
    protein={
        "conf": "path/to/protein.{pdb;gro}",
        "ff":{
            "code": <external_force_field_code>
        }
    },
    custom_ff_path='parent/directory/of/custom.ff'
    ...   
)

by top

Be careful

It is advised to build a single topology file without any include statements. If you want to use those include statements, they MUST BE absolute paths to their corresponded files.

calculate(
    ...
    protein={
        "conf": "path/to/protein.gro",
        "top": "path/to/protein.top",
    },
    ...   
)

membrane

In this case a topology must be generated and provided. This topology can also be obtained from CHARMM-GUI.

Be careful

It is advised to build a single topology file without any include statements. If you want to use those include statements, they MUST BE absolute paths to their corresponded files.

You must always past the PDB despite passing a topology, GRO files are not accepted at the moment.

Getting the membrane.pdb file

For a membrane systems you first need to embed the protein into the membrane. This can easily done with CHARMM-GUI. This is a non-exhaustive list of steps for this process:

Processing on CHARMM-GUI:

• ACE and CT3 terminus

• pH=7

• Run PPM 2.0

• Hexagonal box

• Water thickness: 15

• Length of X and Y: 90 (initial guess)

• Only POPC for a simple membrane or any other lipid composition you would like to use

• Do not include ions

• force field:

  • AMBER:

    • Protein:

    • Lipid:

    • Water: TIP3P

• Input Generation Options: GROMACS

Open gromac/step5_input.gro. This file has the crystal information which is important during the solvation step in BindFLow. Convert to PDB with gmx editconf -f step5_input.gro -o temporal.pdb . The last PDB will also have the crystal info. Split the PDB in POPC and protein.

In PyMOL:

• select popc, resn POPC

• select prot, (polymer.protein or resn ACE or resn NME)

👀 Depending on how you are processing the protein, you might need to change the entry ATOM to HETATM for the CAP groups.

You may also need to manually split the topology into protein and membrane.

calculate(
    ...
    membrane={
        "conf": "path/to/membrane.pdb",
        "top": "path/to/membrane.top",
    }
    ...
)

ligands

openff

Any force field from OpenFF can be acceded by setting its name as code (see that the .offxml extension is kept).

If code is not provided, the default force field for type = "openff" is openff_unconstrained-2.0.0.offxml.

ligand_files = [
    "path/to/ligand1.{mol;sdf}",
    "path/to/ligand2.{mol;sdf}",
    "path/to/ligand3.{mol;sdf}",
    ...
]

ligands = []
for ligand_file in ligand_files:
    ligands.append({
        "conf": ligand_file,
        "ff":{
            "type": "openff",
            "code": "openff_unconstrained-2.0.0.offxml"
        }
    })

calculate(
    ...
    ligands=ligands,
    ...
)

espaloma

It is recommended to use espaloma >= 0.3.1.

If code is not provided, the default force field for type = "espaloma" is espaloma-0.3.1.

ligand_files = [
    "path/to/ligand1.{mol;sdf}",
    "path/to/ligand2.{mol;sdf}",
    "path/to/ligand3.{mol;sdf}",
    ...
]

ligands = []
for ligand_file in ligand_files:
    ligands.append({
        "conf": ligand_file,
        "ff":{
            "type": "espaloma",
            "code": "espaloma-0.3.1"
        }
    })

calculate(
    ...
    ligands=ligands,
    ...
)

gaff

If code is not provided, the default force field for type = "gaff" is gaff-2.11.

ligand_files = [
    "path/to/ligand1.{mol;sdf}",
    "path/to/ligand2.{mol;sdf}",
    "path/to/ligand3.{mol;sdf}",
    ...
]

ligands = []
for ligand_file in ligand_files:
    ligands.append({
        "conf": ligand_file,
        "ff":{
            "type": "gaff",
            "code": "gaff-2.11"
        }
    })

calculate(
    ...
    ligands=ligands,
    ...
)

custom force field

Be careful

It is advised to build a single topology file without any include statements. If you want to use those include statements, they MUST BE absolute paths to their corresponded files.

ligands = [
    {
        "conf": "path/to/ligand1.gro",
        "top": "path/to/ligand1.top"
    },
    {
        "conf": "path/to/ligand2.gro",
        "top": "path/to/ligand2.top"
    },
    {
        "conf": "path/to/ligand3.gro",
        "top": "path/to/ligand3.top"
    },
    ...
]

calculate(
    ...
    ligands=ligands,    
    ...
)

cofactor

openff

Any force field from OpenFF can be acceded by setting its name as code (see that the .offxml extension is kept).

If code is not provided, the default force field for type = "openff" is openff_unconstrained-2.0.0.offxml.

calculate(
    ...
    cofactor={
        "conf": "path/to/cofactor.{mol;sdf}",
        "ff":{
            "type": "openff",
            "code": "openff_unconstrained-2.0.0.offxml"
        }
    },
    ...
)

espaloma

It is recommended to use espaloma >= 0.3.1.

If code is not provided, the default force field for type = "espaloma" is espaloma-0.3.1.

calculate(
    ...
    cofactor={
        "conf": "path/to/cofactor.{mol;sdf}",
        "ff":{
            "type": "espaloma",
            "code": "espaloma-0.3.1"
        }
    },
    ...
)

gaff

If code is not provided, the default force field for type = "gaff" is gaff-2.11.

calculate(
    ...
    cofactor={
        "conf": "path/to/cofactor.{mol;sdf}",
        "ff":{
            "type": "gaff",
            "code": "gaff-2.11"
        }
    },
    ...
)

custom force field

non-water

Be careful

It is advised to build a single topology file without any include statements. If you want to use those include statements, they MUST BE absolute paths to their corresponded files.

calculate(
    ...
    cofactor={
        "conf": "path/to/cofactor.gro",
        "top": "path/to/cofactor.top",
    },    
    ...
)

water-like

Be careful

It is advised to build a single topology file without any include statements. If you want to use those include statements, they MUST BE absolute paths to their corresponded files.

In the case that the cofactor(s) is (are) water-like molecule(s), this should be specified by the keyword is_water = True. In this case, a special treatment is done in BindFlow internally. Here, its settles section (if any) will be changed to TIP3P-like triangular constraints. Check the discussion How to treat specific water molecules as ligand?. Note that this is only possible for TIP3P-like water molecules at the moment.

calculate(
    ...
    cofactor={
        "conf": "path/to/cofactor.gro",
        "top": "path/to/cofactor.top",
        "is_water": True,
    },    
    ...
)

Water models

BindFlow comes with (at present 07.2024) all water models distributed with GROMACS. They are set by the keyword: water_model. E.g.:

calculate(
    ...
    water_model="amber/tip3p"
    ...
)

The structure of the string is force_field_family/water_model; amber/tip3p is the default.

All configuration and topology files are sourced from GROMACS force fields, available at GROMACS GitLab - share/top. These files contain topologies and configurations for water models and ions within three force field families: AMBER, CHARMM, and OPLS-AA.

Note

It is assumed that inside the same family, the non-bonded interactions (for water models and ions) are the same (epsilon and sigma parameters), which is true for the force fields presented in the GROMACS distribution.

The ffnonbonded.itp for each family was taken from:

• AMBER: amber99sb-ildn

• CHARMM: charmm27

• OPLS-AA: oplsaa

These ffnonbonded.itp files were modified to retain only the [ atomtypes ] section, including only atom types related to the water models and ions. This modification prevents potential conflicts with atom-type definitions from user-provided force fields.

The available force fields and their corresponding configuration files are:

amber:
  spc: spc216.gro
  spce: spc216.gro
  tip3p: spc216.gro
  tip4p: tip4p.gro
  tip4pew: tip4p.gro
  tip5p: tip5p.gro
charmm:
  spc: spc216.gro
  spce: spc216.gro
  tip3p: spc216.gro
  tip4p: tip4p.gro
  tip5p: tip5p.gro
  tips3p: spc216.gro
oplsaa:
  spc: spc216.gro
  spce: spc216.gro
  tip3p: spc216.gro
  tip4p: tip4p.gro
  tip4pew: tip4p.gro
  tip5p: tip5p.gro
  tip5pe: tip5p.gro

MDP options modification based on the force field. AMBER and CHARMM36-like force fields example

Some Molecular Dynamic Parameters (MDP) are usually, rather than interchangeable options, parts of each force field derivation and parametrization. So, we should always use those parameters during our simulations. A typical example are AMBER and CHARMM36-like force fields:

AMBER-like force fields

BindFlow use this parameter by default. So, you do not need to modify them.

constraints: all-bonds
cutoff-scheme: Verlet
vdwtype: cutoff
vdw-modifier: Potential-Shift-Verlet
rlist: 1.2
rvdw: 1.0
coulombtype: PME
rcoulomb: 1.0
DispCorr: EnerPres

CHARMM36-like force fields

In this case you should pass to to all the steps these parameters. BindFlow works with AMBER-like force fields by default. In the GROMACS’ docs says:

constraints: h-bonds
cutoff-scheme: Verlet
vdwtype: cutoff
vdw-modifier: force-switch
rlist: 1.2
rvdw: 1.2
rvdw-switch: 1.0
coulombtype: PME
rcoulomb: 1.2
DispCorr: no

Note that dispersion correction should be applied in the case of lipid monolayers, but not bilayers.

Please also note that the switching distance is a matter of some debate in lipid bilayer simulations, and it is dependent to some extent on the nature of the lipid. Some studies have found that an 0.8-1.0 nm switch is appropriate, others argue 0.8-1.2 nm is best, and yet others stand by 1.0-1.2 nm. The user is cautioned to thoroughly investigate the force field literature for their chosen lipid(s) before beginning a simulation!

Final note

Remember to cite properly the main references if you use any of the force fields in your work.