mlpot (c51b1)
PyCHARMM Element doc/mlpot.info 2.0
Machine Learning Potential (MLpot)
By K. Töpfer, L.I. Vazquez-Salazar, O.T. Unke, M. Meuwly
(kai.toepfer@unibas.ch , m.meuwly@unibas.ch)
The MLpot extension in the PyCHARMM toolkit allows using a custom
machine learning-based (ML) interaction potential based on the Asparagus
program package (https://github.com/MMunibas/Asparagus.git) to apply, e.g., a
PhysNet model potential (J. Chem. Theory Comput. 2019, 15, 3678-3693). This
implementation allows the user to run mixed ML/MM simulations akin to the more
widely known QM/MM simulations.
The custom potential energy function does not only provide the
application of a n-dimensional reactive force field for a selected set
of ML atoms, it can also provide conformationally dependent atom
centered charges for computing electrostatic interactions with the
fixed point charges of the other MM atoms in the simulation. The ML-MM
electrostatic interactions are handled by the Asparagus model potential
using the shifted Coulomb cutoff scheme to include the partial derivative
of the ML atom charges by the ML atom position. Alternatively, fixed point
charges as defined in the psf file can be used for the ML part as
well. Note, that the LJ parameter do not adopt changes in the bonding
structure. The implementation of the PhysNet model in CHARMM using the
PyCHARMM API is intended to be as flexible as possible to use for
other ML potential models.
* Syntax | Syntax of MLpot Class
* Description | Description of MLpot Input Parameters
* Requirements | Requirements to Run MLpot Module
* Capabilities | Description of Capabilities
* Examples | Various Examples
* Frequent Errors | Solution for Frequent Errors
* Notes | Developer notes
Syntax of MLpot Class
pycharmm.MLpot(
ml_model,
ml_Z,
ml_selection,
ml_charge=0,
ml_fq=True,
mlmm_ctofnb=None,
mlmm_ctonnb=None,
mlmm_lambda=None,
**kwargs)
Description of MLpot Input Parameters
Keyword - (Type, Default) Purpose
----------------------------------------------------------------------------
ml_model - (asparagus.interface.model_pycharmm.PyCharmm_Calculator)
PyCHARMM model potential calculator object from the
Asparagus program package, which must be initialized prior
to the activation of the pycharmm.MLpot module.
Z - (list) Respective atomic numbers of the ML atoms
which are consistent with the atomic number label in the
trained ML model.
ml_selection - (pycharmm.SelectAtoms object) SelectAtoms object of the
atoms that are treated by the ML potential.
ml_charge - (int, default=0) Total charge of the ML atom selection.
Note: For ml_fq=True, the psf charges
of the ML atom are all set to zero to neglect
electrostatic interaction by the CHARMM force
field. This could be in conflict for a selection
of ML atoms with total charges different than
zero.
ml_fq - (bool, default=True) Use conformational dependent
fluctuating charge prediction of the ML potential model to
compute the ML-MM electrostatic interaction. If False, the
ML-MM electrostatic interaction is captured by the CHARMM
force field using the atomic charges defined in the psf
file.
mlmm_ctofnb - (float, default=None) Maximum non-bonded interaction range
of the ML-MM atom electrostatic interactions. if None, the
ctofnb value defined in CHARMM is used. Note that ML
potential uses the non-bonded list to include
interaction between ML atoms of the central cell with the
ones from the image cells if periodic boundary
conditions are enabled. That requires cutnb larger than the
cutoff range of the PhysNet model and cutnb should be
larger than cutofnb (ctonnb < ctofnb < cutnb).
mlmm_ctonnb - (float, default=None) Non-bonded interaction range threshold
of the ML-MM atom electrostatic interactions at which the
switch-off function is applied additional to the shifted
Coulomb potential scheme. if None, the ctonnb value defined
in CHARMM is used.
mlmm_lambda - (float, default=None) Scaling factor for the ML-MM
electrostatic interaction. If None, no scaling is applied.
Requirements to Run MLpot Module
The MLpot module requires the working Asparagus program package linked
via to by the PYTHONPATH environment variable. It is available at
https://github.com/MMunibas/Asparagus.git and supports the sampling,
training and application of implemented machine learning-based model
potentials using the PyTorch module.
To carry out computations with a ML potential a PyCHARMM calculator
class object of a trained model potential must be initialized in the
PyCHARMM script and passed to the MLpot function. To avoid errors, the
MLpot module should be called in the PyCHARMM script after the molecular
system is fully set up and the non-bonding parameter setup is defined.
There can only be one MLpot module active in the simulation.
Description of Capabilities
The MLpot module is only callable via PyCHARMM. MLpot is compatible with
standard CHARMM machinery for periodic boundaries, non-bonded cutoffs and
constant-pressure simulations. If fluctuating charges are disabled, MLpot
behaves likes a usual custom energy function defined in usersb and Ewald
summation is supported. MLpot reports its energies in two dedicated energy
terms: the internal ML atom potential energy is accumulated in MLPO, and the
ML-MM electrostatic energy in MLEL. Both are printed on the MLPOT> line of the
energy table, and only when they are non-zero, so a run that does not use
MLpot sees the usual energy table. They may be skipped with "SKIP MLPO MLEL".
Note for users of earlier versions: MLpot previously added these energies to
the user energy term (USER) and the electrostatic term (ELEC) respectively.
Scripts that read the ML contributions out of Euser or ELEC, or that disabled
MLpot with "SKIP USER", must be updated to the MLPO and MLEL terms.
If fluctuating charges are enabled, only zero charge ML atom selection are
supported so far and Ewald summation is applicable. As ML-MM electrostatic
interactions are handled by the MLpot module itself, the ML atom charges in
the psf file are set to zero to avoid double counting. In case of Ewald
summation, the near field evaluation of the ML-MM electrostatic interaction
is done by MLpot, but far range interaction between ML-MM atoms becomes zero
due to the zero charges of ML atoms in the psf file.
Technically, all existing bonds, angles, dihedrals and improper angles
between the ML atoms are deleted and ML atom pairs are included to the
non-bonding exclusion list. CHARMM may also deletes all angles, dihedrals
and improper angles including just even one ML atom.
Tests have shown, that the SHAKE bond constraints such as the hydrogen bond
constraints are still working.
No other CHARMM modules are included into the MLpot implementation.
Thus, all modules affecting partially the potential energy including ML atoms
(e.g. BLOCK module) could not be compatible with MLpot and should be
thoroughly tested before production runs are carried out.
If the simulated system only contains atoms defined in the ML atom
selection, CHARMM cannot handle the energy output information
correctly and shows 'NaN' and '-Infinite' values. The simulation
however seems to run correctly.
A conversion of the CHARMM neighbor list to a ML-ML and ML-MM atom
pair list is performed in a CHARMM module (api_func.F90) using module
defined integer arrays of certain size to avoid allocation each time
the function is called. Currently the max size support a ML atom
selection size of 100 atoms (max_Nml) and pair interactions of up to
100000 (max_Npr). Especially the number of pair interaction might
become a source for errors for very large ML and MM atom numbers.
Examples:
One test script for the simulation of a ammonia in water can be
found in the charmm directory (tool/pycharmm/tests/testMLpot.py)
Frequent Errors:
- Error in PyCHARMM when reading system (e.g. read psf file)
Error:
"""
***** LEVEL -3 WARNING FROM <psf_read_formatted> *****
***** atom type not found for atom
******************************************
"""
Possible Solution:
- Different to native CHARMM, PyCHARMM needs an 'ATOMS' block in the
parameter ('*.par' or '*.prm') files. Otherwise it won't recognize
bonding terms. Basically just copy the respective atom definitions from
the topology file ('*.top' or '*.rtf') to the 'ATOMS' block at the top
of the parameter file.
E.g. (examples/FPhOH_water/toppar/pfoh.par):
"""
* Parameters generated by analogy by
* CHARMM General Force Field (CGenFF) program version 1.0.0
*
ATOMS
MASS 61 CG2R61 12.01100 ! 6-mem aromatic C
MASS 65 CG2R66 12.01100 ! 6-mem aromatic carbon bound to F
MASS 24 HGR62 1.00800 ! nonpolar H, neutral 6-mem planar ...
MASS 23 HGR61 1.00800 ! aromatic H
[...]
BONDS
CG2R61 CG2R61 305.00 1.3750 ! PROT benzene, JES 8/25/89
CG2R61 CG2R66 305.00 1.3700 ! NAMODEL difluorotoluene
[...]
"""
- MLpot energy (PhysNet) shows 'NaN'.
Errors:
- None, just 'NaN' in the energy output block for 'MLPO' (in versions
before the dedicated MLpot terms, this appeared under 'Euser')
Possible Solution:
- Check if the atomic number list passed to the MLpot module matches with
the atomic numbers in the training data. A mismatch can likely appear
when using ASE to read a .pdb file and returning the atomic numbers
(e.g. 'Z = atoms.get_atomic_numbers()', see line 102 in
'examples/Criegee/Script_run.py'). ASE interpret the atom name as atom
type and names such as 'CA...' for an alpha-carbon can be seen as atom
type calcium 'Ca'.
- ML treated molecule breaks imediatly or quickly at the start of MD simulation
or even unconstrained minimizations.
Errors:
- ML molecule adopting chemically wrong structure or splits into atoms
Possible Solution:
- In MD simulations, conformations are adopted which are outside the
training samples.
This artifact mostly appears when the ML treated molecule adopts
structures which were not included in the training set. This often
happens for, e.g., ML-MM simulation of ML molecules in MM solution
where different conformation were reached in comparison to gas phase
simulations.
- In structure optimizations, forces of initial structures are high and
lead to bad conformations.
It also can happend during the structure optimization (MINI) when the
new guess is outside the trained conformation because of a huge
potential (and forces) difference between classical force field and the
ML model potential. Normally, large changes or even jumps between
optimization steps are prevented by method and limited by the keyword
'step real' (default: step 0.02). This was not observed in the case of a
'MINI SD' run, that lead to large jumps in bond lengths (up to 0.5 Ang)
after the first step.
Fix: - lower step parameter, e.g., 'MINI SD STEP 0.001'
Machine Learning Potential (MLpot)
By K. Töpfer, L.I. Vazquez-Salazar, O.T. Unke, M. Meuwly
(kai.toepfer@unibas.ch , m.meuwly@unibas.ch)
The MLpot extension in the PyCHARMM toolkit allows using a custom
machine learning-based (ML) interaction potential based on the Asparagus
program package (https://github.com/MMunibas/Asparagus.git) to apply, e.g., a
PhysNet model potential (J. Chem. Theory Comput. 2019, 15, 3678-3693). This
implementation allows the user to run mixed ML/MM simulations akin to the more
widely known QM/MM simulations.
The custom potential energy function does not only provide the
application of a n-dimensional reactive force field for a selected set
of ML atoms, it can also provide conformationally dependent atom
centered charges for computing electrostatic interactions with the
fixed point charges of the other MM atoms in the simulation. The ML-MM
electrostatic interactions are handled by the Asparagus model potential
using the shifted Coulomb cutoff scheme to include the partial derivative
of the ML atom charges by the ML atom position. Alternatively, fixed point
charges as defined in the psf file can be used for the ML part as
well. Note, that the LJ parameter do not adopt changes in the bonding
structure. The implementation of the PhysNet model in CHARMM using the
PyCHARMM API is intended to be as flexible as possible to use for
other ML potential models.
* Syntax | Syntax of MLpot Class
* Description | Description of MLpot Input Parameters
* Requirements | Requirements to Run MLpot Module
* Capabilities | Description of Capabilities
* Examples | Various Examples
* Frequent Errors | Solution for Frequent Errors
* Notes | Developer notes
Syntax of MLpot Class
pycharmm.MLpot(
ml_model,
ml_Z,
ml_selection,
ml_charge=0,
ml_fq=True,
mlmm_ctofnb=None,
mlmm_ctonnb=None,
mlmm_lambda=None,
**kwargs)
Description of MLpot Input Parameters
Keyword - (Type, Default) Purpose
----------------------------------------------------------------------------
ml_model - (asparagus.interface.model_pycharmm.PyCharmm_Calculator)
PyCHARMM model potential calculator object from the
Asparagus program package, which must be initialized prior
to the activation of the pycharmm.MLpot module.
Z - (list) Respective atomic numbers of the ML atoms
which are consistent with the atomic number label in the
trained ML model.
ml_selection - (pycharmm.SelectAtoms object) SelectAtoms object of the
atoms that are treated by the ML potential.
ml_charge - (int, default=0) Total charge of the ML atom selection.
Note: For ml_fq=True, the psf charges
of the ML atom are all set to zero to neglect
electrostatic interaction by the CHARMM force
field. This could be in conflict for a selection
of ML atoms with total charges different than
zero.
ml_fq - (bool, default=True) Use conformational dependent
fluctuating charge prediction of the ML potential model to
compute the ML-MM electrostatic interaction. If False, the
ML-MM electrostatic interaction is captured by the CHARMM
force field using the atomic charges defined in the psf
file.
mlmm_ctofnb - (float, default=None) Maximum non-bonded interaction range
of the ML-MM atom electrostatic interactions. if None, the
ctofnb value defined in CHARMM is used. Note that ML
potential uses the non-bonded list to include
interaction between ML atoms of the central cell with the
ones from the image cells if periodic boundary
conditions are enabled. That requires cutnb larger than the
cutoff range of the PhysNet model and cutnb should be
larger than cutofnb (ctonnb < ctofnb < cutnb).
mlmm_ctonnb - (float, default=None) Non-bonded interaction range threshold
of the ML-MM atom electrostatic interactions at which the
switch-off function is applied additional to the shifted
Coulomb potential scheme. if None, the ctonnb value defined
in CHARMM is used.
mlmm_lambda - (float, default=None) Scaling factor for the ML-MM
electrostatic interaction. If None, no scaling is applied.
Requirements to Run MLpot Module
The MLpot module requires the working Asparagus program package linked
via to by the PYTHONPATH environment variable. It is available at
https://github.com/MMunibas/Asparagus.git and supports the sampling,
training and application of implemented machine learning-based model
potentials using the PyTorch module.
To carry out computations with a ML potential a PyCHARMM calculator
class object of a trained model potential must be initialized in the
PyCHARMM script and passed to the MLpot function. To avoid errors, the
MLpot module should be called in the PyCHARMM script after the molecular
system is fully set up and the non-bonding parameter setup is defined.
There can only be one MLpot module active in the simulation.
Description of Capabilities
The MLpot module is only callable via PyCHARMM. MLpot is compatible with
standard CHARMM machinery for periodic boundaries, non-bonded cutoffs and
constant-pressure simulations. If fluctuating charges are disabled, MLpot
behaves likes a usual custom energy function defined in usersb and Ewald
summation is supported. MLpot reports its energies in two dedicated energy
terms: the internal ML atom potential energy is accumulated in MLPO, and the
ML-MM electrostatic energy in MLEL. Both are printed on the MLPOT> line of the
energy table, and only when they are non-zero, so a run that does not use
MLpot sees the usual energy table. They may be skipped with "SKIP MLPO MLEL".
Note for users of earlier versions: MLpot previously added these energies to
the user energy term (USER) and the electrostatic term (ELEC) respectively.
Scripts that read the ML contributions out of Euser or ELEC, or that disabled
MLpot with "SKIP USER", must be updated to the MLPO and MLEL terms.
If fluctuating charges are enabled, only zero charge ML atom selection are
supported so far and Ewald summation is applicable. As ML-MM electrostatic
interactions are handled by the MLpot module itself, the ML atom charges in
the psf file are set to zero to avoid double counting. In case of Ewald
summation, the near field evaluation of the ML-MM electrostatic interaction
is done by MLpot, but far range interaction between ML-MM atoms becomes zero
due to the zero charges of ML atoms in the psf file.
Technically, all existing bonds, angles, dihedrals and improper angles
between the ML atoms are deleted and ML atom pairs are included to the
non-bonding exclusion list. CHARMM may also deletes all angles, dihedrals
and improper angles including just even one ML atom.
Tests have shown, that the SHAKE bond constraints such as the hydrogen bond
constraints are still working.
No other CHARMM modules are included into the MLpot implementation.
Thus, all modules affecting partially the potential energy including ML atoms
(e.g. BLOCK module) could not be compatible with MLpot and should be
thoroughly tested before production runs are carried out.
If the simulated system only contains atoms defined in the ML atom
selection, CHARMM cannot handle the energy output information
correctly and shows 'NaN' and '-Infinite' values. The simulation
however seems to run correctly.
A conversion of the CHARMM neighbor list to a ML-ML and ML-MM atom
pair list is performed in a CHARMM module (api_func.F90) using module
defined integer arrays of certain size to avoid allocation each time
the function is called. Currently the max size support a ML atom
selection size of 100 atoms (max_Nml) and pair interactions of up to
100000 (max_Npr). Especially the number of pair interaction might
become a source for errors for very large ML and MM atom numbers.
Examples:
One test script for the simulation of a ammonia in water can be
found in the charmm directory (tool/pycharmm/tests/testMLpot.py)
Frequent Errors:
- Error in PyCHARMM when reading system (e.g. read psf file)
Error:
"""
***** LEVEL -3 WARNING FROM <psf_read_formatted> *****
***** atom type not found for atom
******************************************
"""
Possible Solution:
- Different to native CHARMM, PyCHARMM needs an 'ATOMS' block in the
parameter ('*.par' or '*.prm') files. Otherwise it won't recognize
bonding terms. Basically just copy the respective atom definitions from
the topology file ('*.top' or '*.rtf') to the 'ATOMS' block at the top
of the parameter file.
E.g. (examples/FPhOH_water/toppar/pfoh.par):
"""
* Parameters generated by analogy by
* CHARMM General Force Field (CGenFF) program version 1.0.0
*
ATOMS
MASS 61 CG2R61 12.01100 ! 6-mem aromatic C
MASS 65 CG2R66 12.01100 ! 6-mem aromatic carbon bound to F
MASS 24 HGR62 1.00800 ! nonpolar H, neutral 6-mem planar ...
MASS 23 HGR61 1.00800 ! aromatic H
[...]
BONDS
CG2R61 CG2R61 305.00 1.3750 ! PROT benzene, JES 8/25/89
CG2R61 CG2R66 305.00 1.3700 ! NAMODEL difluorotoluene
[...]
"""
- MLpot energy (PhysNet) shows 'NaN'.
Errors:
- None, just 'NaN' in the energy output block for 'MLPO' (in versions
before the dedicated MLpot terms, this appeared under 'Euser')
Possible Solution:
- Check if the atomic number list passed to the MLpot module matches with
the atomic numbers in the training data. A mismatch can likely appear
when using ASE to read a .pdb file and returning the atomic numbers
(e.g. 'Z = atoms.get_atomic_numbers()', see line 102 in
'examples/Criegee/Script_run.py'). ASE interpret the atom name as atom
type and names such as 'CA...' for an alpha-carbon can be seen as atom
type calcium 'Ca'.
- ML treated molecule breaks imediatly or quickly at the start of MD simulation
or even unconstrained minimizations.
Errors:
- ML molecule adopting chemically wrong structure or splits into atoms
Possible Solution:
- In MD simulations, conformations are adopted which are outside the
training samples.
This artifact mostly appears when the ML treated molecule adopts
structures which were not included in the training set. This often
happens for, e.g., ML-MM simulation of ML molecules in MM solution
where different conformation were reached in comparison to gas phase
simulations.
- In structure optimizations, forces of initial structures are high and
lead to bad conformations.
It also can happend during the structure optimization (MINI) when the
new guess is outside the trained conformation because of a huge
potential (and forces) difference between classical force field and the
ML model potential. Normally, large changes or even jumps between
optimization steps are prevented by method and limited by the keyword
'step real' (default: step 0.02). This was not observed in the case of a
'MINI SD' run, that lead to large jumps in bond lengths (up to 0.5 Ang)
after the first step.
Fix: - lower step parameter, e.g., 'MINI SD STEP 0.001'