diff --git a/README.md b/README.md index 09ce92a12..5406747fb 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,27 @@ # DMFF -Differentiable Molecular Force Field + +**DMFF** (**D**ifferentiable **M**olecular **F**orce **F**ield) is a Jax-based python package that provides a full differentiable implementation of molecular force field models. This project aims to establish an extensible codebase to minimize the efforts in force field parameterization, and to ease the force and virial tensor evaluations for advanced complicated potentials (e.g., polarizable models with geometry-dependent atomic parameters). Currently, this project mainly focuses on the molecular systems such as: water, biological macromolecules (peptides, proteins, nucleic acids), organic polymers, and small organic molecules (organic electrolyte, drug-like molecules) etc. We support both the conventional point charge models (OPLS and AMBER like) and multipolar polarizable models (AMOEBA and MPID like). The entire project is backed by the XLA technique in JAX, thus can be "jitted" and run in GPU devices much more efficiently compared to normal python codes. + +The behavior of organic molecular systems (e.g., protein folding, polymer structure, etc.) is often determined by a complex effect of many different types of interactions. The existing organic molecular force fields are mainly empirically fitted and their performance relies heavily on error cancellation. Therefore, the transferability and the prediction power of these force fields are insufficient. For new molecules, the parameter fitting process requires essential manual intervention and can be quite cumbersome. In order to automate the parametrization process and increase the robustness of the model, it is necessary to apply modern AI techniques in conventional force field development. This project serves for this purpose by utilizing the automatic differentiable programming technique to develop a codebase, which allows a more convenient incorporation of modern AI optimization techniques. It also helps the realization of many exciting functions including (but not limited to): hybrid machine learning/force field models and parameter optimization based on trajectory. + +## User Guide + ++ [1. Introduction](user_guide/introduction.md) ++ [2. Installation](user_guide/installation.md) ++ [3. Compute energy and forces](user_guide/compute.md) ++ [4. Compute gradients with auto differentiable framework](user_guide/auto_diff.md) ++ [5. Theories](user_guide/theory.md) ++ [6. Introduction to force field xml files](user_guide/xml_spec.md) + +## Developer Guide ++ [1. Introduction](dev_guide/introduction.md) ++ [2. Architecture](dev_guide/arch.md) ++ [3. Convention](dev_guide/convention.md) + +## Modules ++ [1. ADMP](modules/admp.md) + + +## Support and Contribution + +Please visit our repository on [GitHub](https://github.com/deepmodeling/DMFF) for the library source code. Any issues or bugs may be reported at our issue tracker. All contributions to DMFF are welcomed via pull requests! diff --git a/docs/about.md b/docs/about.md index e69de29bb..835deeb6f 100644 --- a/docs/about.md +++ b/docs/about.md @@ -0,0 +1,3 @@ +# Lisense + +# Contributor \ No newline at end of file diff --git a/docs/admp/README.md b/docs/admp/README.md deleted file mode 100644 index 9ca538f1a..000000000 --- a/docs/admp/README.md +++ /dev/null @@ -1,68 +0,0 @@ -# ADMP - -Automatic Differentiable Multipolar Polarizable (ADMP) force field calculator. - -This module provides an auto-differentiable implementation of multipolar polarizable force fields, that resembles the behavior of [MPID](https://github.com/andysim/MPIDOpenMMPlugin) plugin of OpenMM. Supposedly, this module is developed for the following purposes: - -1. Achieving an easy calculation of force and virial of the multipolar polarizable forcefield. -2. Allowing fluctuating (geometric-dependent) multipoles/polarizabilities in multipolar polarizable potentials. -3. Allowing the calculation of derivatives of various force field parameters, thus achieving a more systematic and automatic parameter optimization scheme. - -The module is based on [JAX](https://github.com/google/jax) and [JAX-MD](https://github.com/google/jax-md) projects. - - - -## Installation - -### Dependencies - -ADMP module depends on the following packages, install them before using ADMP: - -1. Install [jax](https://github.com/google/jax) (pick the correct cuda version, see more details on their installation guide): - - ```bash - pip install jax[cuda11_cudnn82] -f https://storage.googleapis.com/jax-releases/jax_releases.html - ``` - -2. Install [jax-md](https://github.com/google/jax-md) : - - ```bash - pip install jax-md --upgrade - ``` - - ADMP currently relies on the space and partition modules to provide neighbor list - -3. Install ADMP: - - ADMP is a pure python module, just simply put it in your $PYTHONPATH. - - ```bash - export PYTHONPATH=$PYTHONPATH:/path/to/admp - ``` - - - -## Settings - -In `admp/settings.py`, you can modify some global settings, including: - -**PRECISION**: single or double precision - -**DO_JIT**: whether do jit or not. - - - -## Example - -We provide a MPID 1024 water box example. In water_1024 and water_pol_1024, we show both the nonpolarizable and the polarizable cases. - -```bash -cd ./examples/water_1024 -./run_admp.py - -cd ./examples/water_pol_1024 -./run_admp.py -``` - -if `DO_JIT = True`, then the first run would be a bit slow, since it tries to do the jit compilation. Further executions of `get_forces` or `get_energy` should be much faster. - diff --git a/docs/admp/frontend.md b/docs/admp/frontend.md new file mode 100644 index 000000000..fa456fb73 --- /dev/null +++ b/docs/admp/frontend.md @@ -0,0 +1,85 @@ +# Implemented Short-Range Potentials in DMFF + +DMFF already provides the frontends for many frequently-used short-range potentials. +We document them in this page: + +## SlaterExForce + +Slater-type short-range repulsive force +Ref: jctc 12 3851 + +Example: +```xml + + + + + + + +``` + +Formula: + +$$ +\displaylines{ +E = \sum_{ij} {A_{ij}P(B_{ij}, r)\exp(-B_{ij}r)} \\ +P(B_{ij}, r) = \frac{1}{3}B_{ij}^2 r^2 + B_{ij} r + 1 \\ +A_{ij} = A_i A_j \\ +B_{ij} = \sqrt{B_i B_j} +} +$$ + +## QqTtDampingForce + +Charge-Charge Tang-Tonnies damping force, used in combination with normal PME + +Example: +```xml + + + + + + + +``` + +Formula: + +$$ +\displaylines{ +E = \sum_{ij} {- e^{-B_{ij} r} (1 + B_{ij}r) \frac{q_i q_j}{r}} \\ +B_{ij} = \sqrt{B_i B_j} +} +$$ + +## SlaterDampingForce + +Slater-type damping function for dispersion. Used in combination with the normal dispersion PME +Ref: jctc 12 3851 + +Example: +```xml + + + + + + + +``` + +Formula: + +$$ +\displaylines{ +E = -\sum_{n=6,8,10} {\sum_{ij} {f_n(x)\frac{\sqrt{C_i^n C_j^n}}{r^n}}} \\ +f_n(x) = - e^{-B_{ij} x}\sum_{k=0}^n {\frac{(B_{ij} x)^k}{k!}} \\ +x = B_{ij}r - \frac{2(B_{ij}r)^2 + 3B_{ij}r}{(B_{ij}r)^2 + 3B_{ij}r + 3} \\ +B_{ij} = \sqrt{B_i B_j} +} +$$ diff --git a/docs/admp/readme.md b/docs/admp/readme.md new file mode 100644 index 000000000..b0e157372 --- /dev/null +++ b/docs/admp/readme.md @@ -0,0 +1,258 @@ +# ADMP + +Automatic Differentiable Multipolar Polarizable (ADMP) force field calculator. + +This module provides an auto-differentiable implementation of multipolar polarizable force fields, that resembles the behavior of [MPID](https://github.com/andysim/MPIDOpenMMPlugin) plugin of OpenMM. Supposedly, this module is developed for the following purposes: + +1. Achieving an easy calculation of force and virial of the multipolar polarizable forcefield. +2. Allowing fluctuating (geometric-dependent) multipoles/polarizabilities in multipolar polarizable potentials. +3. Allowing the calculation of derivatives of various force field parameters, thus achieving a more systematic and automatic parameter optimization scheme. + +The module is based on [JAX](https://github.com/google/jax) and [JAX-MD](https://github.com/google/jax-md) projects. + + +## Settings + +In `admp/settings.py`, you can modify some global settings, including: + +**PRECISION**: single or double precision + +**DO_JIT**: whether do jit or not. + +In admp/pme.py, you can also modify the `DEFAULT_THOLE_WIDTH` variable +(You can directly set `dmff.admp.pme.DEFAULT_THOLE_WIDTH` in your code) + + +## Example + +We provide a polarizable 1024 water box example in the water_fullpol folder: + +```bash +cd ./examples/water_fullpol +./run.py +``` + +if `DO_JIT = True`, then the first run would be a bit slow, since it tries to do the jit compilation. Further executions of `get_forces` or `get_energy` should be much faster. + +Following this example, we will introduce the frontend and the backend of the ADMP module. + +## Frontend + +In the `forcefield.xml` file, you can find the frontends of two force field components: + +### ADMPDispForce + +The corresponding XML node is like: + +```xml + + + + +``` + +It computes the following function: + +$$ +\displaylines{ +E = \sum_{i + + + + + + +``` + +In here, the `mScales`, `pScales`, and `dScales` are explained in the [theory](../user_guide/theory.md) page. `lmax` specifies +the highest order of multipoles (`lmax=2` means up to quadrupole). All input parameters that follows are in `nm` and `kJ/mol`. +Currently, we only support the isotropic polarizabilities, meaning the `XX`, `YY`, and `ZZ` components will be averaged. + +The multipole moments are specified in local frames with Cartesian representations. One needs to be careful about the convention +of the Cartesian tensor: in the frontend of ADMP, we adopt the openmm convention, the quadrupole value of which is 3 times smaller +than the convention adopted in Anthony's book. + +In `run.py`, when creating the potential function, several key parameters are noted: +```python +potentials = H.createPotential(pdb.topology, nonbondedCutoff=rc*angstrom, nonbondedMethod=CutoffPeriodic, ethresh=1e-4) +``` +* `ethresh`: this is the energy threshold for PME +* `rc`: the cutoff distance. `rc`, `ethresh`, and `box` together, determine the $K_{max}$ and $\kappa$ + (please see [theory](../user_guide/theory.md)). + Note that the `rc` variable in here is only used to determine the PME settings. The user has to make sure the `rc` + value used in here is the same as the one used in neighbor list construction. +* `nonbondedMethod`: Currently two methods are supported: `CutoffPeriodic` and `PME` (default). When `CutoffPeriodic` + is used, PME is turned off by setting $\kappa=0$ and removing the reciprocal space contribution. + + +## Backend + +The backend of ADMP can be invoked independently without the frontend API, which is much more flexible. +It can be utilized to implement fancy force fields with fluctuating atomic parameters. +***Note that all backend functions are assuming `angstrom` unit, instead of `nm`.*** This is different to the frontend! + +The examples for backend are: examples/water_1024 and examples/water_pol_1024 + +There are mainly three calculators implemented in ADMP backend: + +### ADMPPmeForce + +The `ADMPPmeForce` is initialized as: + +```python +pme_force = ADMPPmeForce(box, axis_type, axis_indices, covalent_map, rc, ethresh, lmax, lpol=True, lpme=True, steps_pol=None) +``` + +The inputs are: + +* `box`: the 3*3 box matrix. Vectors are arranged in rows, in angstrom. + +* `axis_type`: the axis type for the local frame definition for each atom: +```python +ZThenX = 0 +Bisector = 1 +ZBisect = 2 +ThreeFold = 3 +ZOnly = 4 +NoAxisType = 5 +LastAxisTypeIndex = 6 +``` + +* `axis_indices`: indices for local axis definitions, see introduction to local frame in [theory](../user_guide/theory.md) + +* `covalent_map`: a $N\times N$ ($N$ being the number of atoms) matrix of covalent spacings between atoms. + $n$ means the two atoms are $n$ bonds away, and 0 means the two atoms are considered as not bonded topologically. + +* `rc`: cutoff distance in real space. Only used to determine the PME settings. + +* `ethresh`: energy threshold for PME. + +* `lmax`: max L for multipoles. + +* `lpol`: whether turn on polarization? + +* `lpme`: wether turn on PME? + +* `steps_pol`: specifies number of SCF steps when solving induced dipoles. If set to None, then the SCF iteration will stop + when the field is less than `dmff.admp.settings.POL_CONV`. Otherwise, exact `steps_pol` number of iterations will be performed. + ***This tag is important, and needs to be set if you want to jit the function externally.*** + +Once the `pme_force` is initialized, the differentiable calculator can be called as: + +```python +# if lpol = False +E = pme_force.get_energy(positions, box, pairs, Q_local, pol, tholes, mScales, pScales, dScales, U_init=U) +# if lpol = True +E = pme_force.get_energy(positions, box, pairs, Q_local, mScales) +``` + +Important inputs include: + +* `Q_local`: spherical harmonic multipole moments in local frame (in Angstrom). + +* `pol`: polarizability for each atom + +* `tholes`: thole width for each atom + +* `mScales`, `pScales`, `dScales`: topological scaling factors + +* `U_init`: initial values for induced dipoles. + + +### ADMPDispPmeForce + +This force computes the long range dispersion interactions with C6, C8, and C10 terms. +It is initialized and called as: + +```python +disp_force = ADMPDispPmeForce(box, covalent_map, rc, ethresh, pmax, lpme=True) +E = disp_force.get_energy(positions, box, pairs, c_list, mScales) +``` + +Most parameters are similar to the PME force, several other important parameters include: + +* `pmax`: maximum power of dispersion, usually 10. + +* `c_list`: the ***square root*** of dispersion coefficients. (N*3) array, that is, + sqrt(C6), sqrt(C8), and sqrt(C10) of each atom. + + +### Short-Range Pairwise Force + +DMFF provides a simple approach to raise a distance-dependent pair interaction kernel into +a many-body potential function. Suppose you want to implement a function looks like this: + +$$ +E = \sum_{ij} {f(r_{ij}; a_i, a_j, b_i, b_j, c_i, c_j)} +$$ + +Where $a,b,c$ are atomic parameters. + +Then you need to define a pairwise kernel looks like: + +```python +from jax import vmap, jit + +def f_kernel(r, m, a_i, a_j, b_i, b_j, c_i, c_j): + # do calculations, get E + return E * m + +f_kernel = vmap(jit(f_kernel)) +``` + +Then raise it to a many-body potential for a particular system: + +```python +potential_fn = generate_pairwise_interaction(f_kernel, covalent_map) +# a_list, b_list, c_list are atomic parameter lists +E = potential_fn(positions, box, pairs, mScales, a_list, b_list, c_list) +``` + +The resulted `potential_fn` function will take care the following business for you: + +1. Topological scaling based on `covalent_map` and `mScales`; + +2. Assign parameters for each interacting pair; + +We can use this to immplement short-range repulsion and short-range damping quite easily. + + + + + diff --git a/docs/assets/DMFF_arch.md b/docs/assets/DMFF_arch.md new file mode 100644 index 000000000..e0bd17511 --- /dev/null +++ b/docs/assets/DMFF_arch.md @@ -0,0 +1,73 @@ +```mermaid +flowchart TB +E(["atomic & topological params"]) +subgraph G1["Parser & Typification"] +A["force field xml"] --> B(["parseElement"]) +B --> C[Generators with forcefield params] +C --> D(["createPotential"]) +D --> E +end +subgraph G2["ADMP Calculators"] +E --> |init|F(["General Pairwise Calculator"]) +E --> |init|G(["Multipole PME Caculator"]) +E --> |init|H(["Dispersion PME Calculator"]) +end +subgraph G3["Classical Calculators"] +E --> |init|K["Intramol Calculators"] +E --> |init|L["Intermol Calculators"] +end +J["pairs (neighbor list from jax-md)"] --> |input|I +F --> |return|I(["potential(pos, box, pairs, params)"]) +G --> |return|I +H --> |return|I +K --> |return|I +L --> |return|I +C --> |differentiable generator.params|I +``` + + + +```mermaid +graph +A("kernel(dr, m, p0i, p0j, p1i, p1j, ...)") --> D([generate_pairwise_interaction]) +B[covalent_map] --> D +C[static variables] --> D +D --> E(["pair_energy(pos, box, pairs, p0list, p1list, ...)"]) +``` + +```mermaid +classDiagram + class ADMPPmeForce + ADMPPmeForce : +axis_type + ADMPPmeForce : +axis_indices + ADMPPmeForce : +rc + ADMPPmeForce : +ethresh + ADMPPmeForce : +kappa, K1, K2, K3 + ADMPPmeForce : +pme_order + ADMPPmeForce : +covalent_map + ADMPPmeForce : +lpol + ADMPPmeForce : +n_atoms + ADMPPmeForce : +__init__(self, box, axis_type, axis_indices, covalent_map, rc, ethresh, lmax, lpol=False) + ADMPPmeForce : +update_env(attr, val) + ADMPPmeForce : +get_energy(pos, box, pairs, Q_local, pol, tholes, mScales, pScales, dScales) + +``` + + + +``` Mermaid +classDiagram + class ADMPDispPmeForce + ADMPDispPmeForce : +kappa, K1, K2, K3 + ADMPDispPmeForce : +pme_order + ADMPDispPmeForce : +rc + ADMPDispPmeForce : +ethresh + ADMPDispPmeForce : +pmax + ADMPDispPmeForce : +covalent_map + ADMPDispPmeForce : +__init__(self, box, covalent_map, rc, ethresh, pmax) + ADMPDispPmeForce : +update_env(attr, val) + ADMPDispPmeForce : +get_energy(positions, box, pairs, c_list, mScales) +``` + + + diff --git a/docs/assets/arch.png b/docs/assets/arch.png new file mode 100644 index 000000000..4d4a39c08 Binary files /dev/null and b/docs/assets/arch.png differ diff --git a/docs/assets/generator.svg b/docs/assets/generator.svg new file mode 100644 index 000000000..497df01ee --- /dev/null +++ b/docs/assets/generator.svg @@ -0,0 +1 @@ +

Hamiltonian

暴露给用户，继承自OpenMM Forcefield

+ _generators: dict

+ createPotential(topology, **args)

   读入OpenMM Topology，使用

   ResidueTemplate创建SystemData，

   而后构造势函数

Generator

存储力场参数，并负责从拓扑结构生成

整体的势函数

+ param: jax array

   力场参数信息

+ meta: dict

    势函数构建所需信息

+ __init__()

+ parseElement(xml node)

   从xml node中读取参数

+ createForce(SystemData)

   使用拓扑信息创建JAX势函数

+ renderXML()

   使用XML结构化存储力场参数

 \ No newline at end of file diff --git a/docs/assets/opemm_workflow.svg b/docs/assets/opemm_workflow.svg new file mode 100644 index 000000000..17e41aabc --- /dev/null +++ b/docs/assets/opemm_workflow.svg @@ -0,0 +1 @@ +
创建app.forcefield.parser字典
将所有Generator的parseElement方法注册到app.forcefield.parser中
初始化Residue Template部分
调用parser字典中所有的parseElement方法解析XML
将每个调用到的Generator注册到ff object中
用Residue Template对topology进行匹配以分配AtomType
调用所有已注册的Generator的createForce方法来创建对应topology的势函数
import openmm.app as app
ff = app.ForceField("forcefield.xml")
system = ff.createSystem(topology, args)
流程结束
 \ No newline at end of file diff --git a/docs/dev_guide/arch.md b/docs/dev_guide/arch.md new file mode 100644 index 000000000..0204000df --- /dev/null +++ b/docs/dev_guide/arch.md @@ -0,0 +1,352 @@ +# Architecture of DMFF + +![arch](../assets/arch.png) + +The overall architechture of DMFF can be divided into two parts: 1. parser & typing and 2. calculators. +We usually refer to the former as the *frontend* and the latter as the *backend* for ease of description. + +DMFF introduces different forms of force fields in a modular way. For each force field form +(in OpenMM, each form is also called a `Force`), there is a frontend module and a backend module. +The frontend module is responsible for input file parsing, molecular topology construction, atom typification, +and dispatching the forcefield parameters into the atomic parameter; The backend module is the calculation kernel, +which calculates the energy & force of the system, using particle positions and atomic parameters as inputs. + +In the design of the frontend modules, DMFF reuses the frontend parser from OpenMM for topology analysis. +The core class in frontend is the `Generator` class, which should be defined for each force field form. +All frontend `Generators` are currently put in `api.py` and are called by the top-level `Hamiltonian` class. + +The backend module is usually an automatically differentiable calculator built with Jax. + +The structures of the frontend and the backend modules will be introduced in detail in below. + +## How Frontend Works + +Frontend modules are stored in `api.py`. `Hamiltonian` class is the top-level class exposed to users by DMFF. +`Hamiltonian` class reads the path of the XML file, parses the XML file, and calls different frontend modules +according to the `*Force` tags found in XML. The frontend generator has the same form as the OpenMM's generators +[forcefield.py](https://github.com/openmm/openmm/blob/master/wrappers/python/openmm/app/forcefield.py). +The `Generator` class parses the corresponding XML tag to obtain the force field parameters, +then use the parameters to initialize the backend calculator. +It also provides the interface to wrap the backend calculators into potential functions, +which are eventually returned to the users. + +When users use DMFF, the only thing they need to do is to initilize the the `Hamiltonian` class. +In this process, `Hamiltonian` will automatically parse and initialize the corresponding potential function +according to the tags in XML. The call logic is shown in the following chart. +The box represents the command executed in Python script, +and the rounded box represents the internal operation logic of OpenMM when executing the command. + +![openmm_workflow](../assets/opemm_workflow.svg) + +### Hamiltonian Class + +The `Hamiltonian` class is the top-level frontend module, which inherits the +[forcefield class](https://github.com/openmm/openmm/blob/master/wrappers/python/openmm/app/forcefield.py) in OpenMM. +It is responsible for parsing the XML force field files and generating potential functions to calculate system energy +with the given topology information. First, the usage of the `Hamiltonian` class is given: + + +```python +H = Hamiltonian('forcefield.xml') +app.Topology.loadBondDefinitions("residues.xml") +pdb = app.PDBFile("waterbox_31ang.pdb") +rc = 4.0 +# generator stores all force field parameters +generators = H.getGenerators() +disp_generator = generators[0] +pme_generator = generators[1] + +potentials = H.createPotential(pdb.topology, nonbondedCutoff=rc*unit.angstrom) +# pot_fn is the actual energy calculator +pot_disp = potentials[0] +pot_pme = potentials[1] +``` + +`Hamiltonian` class performs the following operations during instantiation: + +* read Residue tag in XML, generate Residue template; +* read AtomTypes tag, store AtomType of each atom; +* for each Force tag, call corresponding `parseElement` method in `app.forcefield.parser` to parse itself, and register `generator`. + +`app.forcefield.parser` is a `dict`, the keys are Force tag names, and the values are the `parseElement` method +of the corresponding `generator`. When `Hamiltonian` parses the XML file, it will use the tag name to look up the +corresponding `parseElement` method, then calls it to initialize the `generator` instance, which stores the raw +parameters from the XML file. You can access all the generators by the `getGenerators()` method in Hamiltonian. + + +### Generator Class + + +The generator class is in charge of input file analysis, molecular topology construction, atom classification, +and expanding force field parameters to atomic parameters. It is a middle layerthat links `Hamiltonian` with the backend. +See the following chart for its design logic: + +![generator](../assets/generator.svg) + +In a custom generator one must define the following methods: + + +* @staticmethod parseElement(element, hamiltonian): OpenMM uses `element.etree` to parse tags in XML file, +and the first argument `element` is an `Element` object. For instance, if there were a section in XML file +that defines a harmonic bond stretching potential: + +```xml + + + + +``` + +This input will be converted into an `Element` object first and be passed to the `HarmonicBondJaxGenerator.parseElement` +method, which should have been registered in `app.forcefield.parsers["HarmonicBondForce"]`. +The developer needs to define the `parseElement` method to parse the input `Element` object and to initialize the generator itself. +For example, in this case, you can use `element.findall("Bond")` to get an iterator that iterates over all the tags. +Then you can use `.attrib` to access all the properties (`type1`, `type2`, `length`, and `k`) within each tag in `dict` format. +These properties, as force field parameters, can be classified into two categories. The first category is differentiable parameters +(such as `length` and `k`), which may be subject to further optimization. By design, these parameters should be fed into the potential +function as explicit input arguments. Therefore, these parameters should be gathered in a `dict` object named `params`, which is then +saved as an attribute of the `generator`. The second category is non-differentiable variables (such as `type1` and `type2`): it is +unlikely that you are going to optimize them, so they are also called *static* variables. These variables will be used by the potential +function implicitly and not be exposed to the users. Therefore, you may save them in `generator` at will, as long as they can be +accessed by the potential function later. + + +* `createForce(self, system, data, nonbondedMethod, *args)` pre-process the force field parameters from XML, and use them to initialize +the backend calculators, then wrap the calculators using a potential function and returns it. +`System` and `data` are given by OpenMM's forcefield class, which store topology/atomType information (for now you need to use +debug tool to access). Bear in mind that we should not break the derivative chain from the XML raw data (force field parameters) +to the per-atom properties (atomic parameters). So we should do the parameter dispatch (i.e., dispatching force field parameters to +each atom) within the returned potential function. Therefore, in `createForce`, we usually only construct the atomtype index map using +the information in `data`, but do not dispatch parameters! The atomtype map will be used in potential function implicitly, to dispatch +parameters. + +Here is an example: + +```python +map_atomtype = np.zeros(n_atoms, dtype=int) + +for i in range(n_atoms): + atype = data.atomType[data.atoms[i]] + map_atomtype[i] = np.where(self.types == atype)[0][0] +``` + +Finally, we need to bind the calculator's compute function to `self._jaxPotential`, which is the final potential function (`potential_fn`) +returned to users: + + +```python + +def potential_fn(positions, box, pairs, params): + return bforce.get_energy( + positions, box, pairs, params["k"], params["length"] + ) + +self._jaxPotential = potential_fn +``` + +The `potential_fn` function only takes `(positions, box, pairs, params)` as explicit input arguments. All these arguments except +`pairs` (neighbor list) should be differentiable. Non differentiable parameters are passed into it by closure (see code convention section). +Meanwhile, if the generator needs to initialize multiple calculators (e.g. `NonBondedJaxGenerator` will call both `LJ` and `PME` calculators), +`potential_fn` should return the summation of the results of all calculators. + +Here is a pseudo-code of the frontend module, demonstrating basic API and method + +```python +from openmm import app + +class SimpleJAXGenerator: + + def __init__(self, hamiltonian): + self.ff = hamiltonian + self.params = None + self._jaxPotential = None + init_other_attributes_if_needed + + @staticmethod + def parseElement(element, hamiltonian): + parse_xml_element + generator = SimpleGenerator(hamiltonian, args_from_xml) + hamiltonian.registerGenerator(generator) + + def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): + generate_constraints_if_needed + # Create JAX energy function from system information + create_jax_potential + self._jaxPotential = jaxPotential + + def getJaxPotential(self, data, **args): + return self._jaxPotential + + def renderXML(self): + render_xml_forcefield_from_params + + +app.parsers["SimpleJAXForce"] = SimpleJAXGenerator.parseElement + +class Hamiltonian(app.ForceField): + + def __init__(self, **args): + super(app.ForceField, self).__init__(self, **args) + self._potentials = [] + + def createPotential(self, topology, **args): + system = self.createSystem(topology, **args) + load_constraints_from_system_if_needed + # create potentials + for generator in self._generators: + potentialImpl = generator.getJaxPotential(data) + self._potentials.append(potentialImpl) + return [p for p in self._potentials] +``` + +And here is a HarmonicBond potential implement: + +```python +class HarmonicBondJaxGenerator: + def __init__(self, hamiltonian): + self.ff = hamiltonian + self.params = {"k": [], "length": []} + self._jaxPotential = None + self.types = [] + + def registerBondType(self, bond): + + types = self.ff._findAtomTypes(bond, 2) + # self.ff._findAtomTypes is a function implemented in OpenMM to patch + # atom types. The first argument is xml element. The second argument is + # the number of types needed to be patched. + # The return of this function is: + # [[atype1, atype2, ...], [atype3, atype4, ...], ...] + # When patching by atom types, the function would return a list with + # patched atom types. When patching by atom classes, the function would + # return a list with all the atom types patched to the class. + + self.types.append(types) + self.params["k"].append(float(bond["k"])) + self.params["length"].append(float(bond["length"])) # length := r0 + + @staticmethod + def parseElement(element, hamiltonian): + # Work with xml tree. Element is the node of forcefield. + # Use element.findall and element.attrib to get the + # children nodes and attributes in the node. + + generator = HarmonicBondJaxGenerator(hamiltonian) + hamiltonian.registerGenerator(generator) + for bondtype in element.findall("Bond"): + generator.registerBondType(bondtype.attrib) + + def createForce(self, system, data, nonbondedMethod, nonbondedCutoff, args): + # jax it! + for k in self.params.keys(): + self.params[k] = jnp.array(self.params[k]) + self.types = np.array(self.types) + + n_bonds = len(data.bonds) + # data is the data structure built by OpenMM, saving topology information of the system. + # The object maintains all the bonds, angles, dihedrals and impropers. + # And it also maintains the atomtype of each particle. + # Use data.atoms, data.bonds, data.angles, data.dihedrals, data.impropers + # to get the atom types. + map_atom1 = np.zeros(n_bonds, dtype=int) + map_atom2 = np.zeros(n_bonds, dtype=int) + map_param = np.zeros(n_bonds, dtype=int) + for i in range(n_bonds): + idx1 = data.bonds[i].atom1 + idx2 = data.bonds[i].atom2 + type1 = data.atomType[data.atoms[idx1]] + type2 = data.atomType[data.atoms[idx2]] + ifFound = False + for ii in range(len(self.types)): + if (type1 in self.types[ii][0] and type2 in self.types[ii][1]) or ( + type1 in self.types[ii][1] and type2 in self.types[ii][0] + ): + map_atom1[i] = idx1 + map_atom2[i] = idx2 + map_param[i] = ii + ifFound = True + break + if not ifFound: + raise BaseException("No parameter for bond %i - %i" % (idx1, idx2)) + + # HarmonicBondJaxForce is the backend class to build potential function + bforce = HarmonicBondJaxForce(map_atom1, map_atom2, map_param) + + # potential_fn is the function to call potential, in which the dict self.params + # is fed to harmonic bond potential function + + def potential_fn(positions, box, pairs, params): + return bforce.get_energy( + positions, box, pairs, params["k"], params["length"] + ) + + self._jaxPotential = potential_fn + # self._top_data = data + + def getJaxPotential(self): + return self._jaxPotential + +# register all parsers +app.forcefield.parsers["HarmonicBondForce"] = HarmonicBondJaxGenerator.parseElement +``` + +## How Backend Works + +### Force Class + +Force class is the backend module that wraps the calculator function. +It does not rely on OpenMM and can be very flexible. For instance, +the Force class of harmonic bond potential is shown below as an example. + +```python +def distance(p1v, p2v): + return jnp.sqrt(jnp.sum(jnp.power(p1v - p2v, 2), axis=1)) + + +class HarmonicBondJaxForce: + def __init__(self, p1idx, p2idx, prmidx): + self.p1idx = p1idx + self.p2idx = p2idx + self.prmidx = prmidx + self.refresh_calculators() + + def generate_get_energy(self): + def get_energy(positions, box, pairs, k, length): + p1 = positions[self.p1idx] + p2 = positions[self.p2idx] + kprm = k[self.prmidx] + b0prm = length[self.prmidx] + dist = distance(p1, p2) + return jnp.sum(0.5 * kprm * jnp.power(dist - b0prm, 2)) + + return get_energy + + def update_env(self, attr, val): + """ + Update the environment of the calculator + """ + setattr(self, attr, val) + self.refresh_calculators() + + def refresh_calculators(self): + """ + refresh the energy and force calculators according to the current environment + """ + self.get_energy = self.generate_get_energy() + self.get_forces = value_and_grad(self.get_energy) +``` + +The design logic for the `Force` class is: it saves the *static* variables inside the class as +the *environment* of the real calculators. Examples of the static environment variables include: +the $\kappa$ and $K_{max}$ in PME calculators, the covalent_map in real-space calculators etc. +For a typical `Force` class, one needs to define the following methods: + +* `__init__`: The initializer, saves all the *static* variables. +* `update_env(self, attr, val)`: updates the saved *static* variables, and refresh the calculators +* `generate_get_energy(self)`: generate a potential calculator named `get_energy` using the current + environment. +* `refresh_calculators(self)`: refresh all calculators if environment is updated. + +In ADMP, all backend calculators only take atomic parameters as input, so they can be invoked +independently in hybrid ML/force field models. The dispatch of force field parameters is done +in the `potential_fn` function defined in the frontend. diff --git a/docs/dev_guide/convention.md b/docs/dev_guide/convention.md new file mode 100644 index 000000000..6aa83f78a --- /dev/null +++ b/docs/dev_guide/convention.md @@ -0,0 +1,22 @@ +# Code Convention + +In this section, you will learn: + - How is DMFF organized + - + +## code organization + +The root directory of DMFF has following sub-directories: + + - `dmff`: source code of project + - `docs`: documents in markdown + - `examples`: examples can be run independently + - `tests`: unit and integration tests + +Under the `dmff`, there are several files and sub-directory: + + - `api.py`: store all the frontend modules + - `settings.py`: global settings + - `utils.py`: helper functions + — each sub-directory represents a set of potential form, e.g. `admp` is Automatic Differentiable Multipolar Polarizable, `classical` is differentiable GAFF forcefield. + diff --git a/docs/dev_guide/introduction.md b/docs/dev_guide/introduction.md new file mode 100644 index 000000000..864be0e1e --- /dev/null +++ b/docs/dev_guide/introduction.md @@ -0,0 +1,21 @@ +# 1. Introduction + +In the developer's guide, you will learn: + ++ Architecture of DMFF ++ Code convention ++ Implement a new function form ++ Write docs ++ Checklist before PR + +DMFF aims to establish a force field calculation and parameter fitting framework based on the automatic differentiation technique. For different forms of force field, DMFF calculates energy and energy gradients through different sub-modules. DMFF ensures sufficient decoupling and modularization at the beginning of the design so new modules can be added easily. + +In the *Architecture of DMFF* section, the overall design of DMFF will be introduced carefully. We will talk about how the parameters are loaded from the XML file and how they are organized for the following energy and gradient calculations. This work is primarily carried out by the `Generator` class. Then we will explain the heavy-lifting `calculators`, which are pure functions that take atomic positions and force field parameters as inputs and compute the energy. + +When developing DMFF, the developers need to obey some programming styles and standards. These are explained in the *Code convention* section. These conventions will help the code reviewers and other developers to understand and maintain the project. + +In the *Implementation of New Potentials* section, we will show how to implement a new calculator for a new potential. Before you turn your equations into code, this section will introduce the specs that your code should follow, making it consistent with other force field sub-modules in DMFF. + +In the *Document Writing* section, we will talk about how to write docs for your new module. DMFF is a collection of force field calculators, and each force field calculator has its own parameters and may be invoked in different ways. Therefore, to make it easy to use and easy to maintain, the developers are required to document the theroy behind the code and the user interface of the code before PR. + +In the *Checklist before PR* section is what you should do before you submit your PR to Github. In this section, you will learn how to write unit tests, check format, and add proper commentsin your code. diff --git a/docs/dev_guide/profile.md b/docs/dev_guide/profile.md new file mode 100644 index 000000000..14a981b77 --- /dev/null +++ b/docs/dev_guide/profile.md @@ -0,0 +1,2 @@ +# How to profile + diff --git a/docs/dev_guide/write_docs.md b/docs/dev_guide/write_docs.md new file mode 100644 index 000000000..3c6019c1d --- /dev/null +++ b/docs/dev_guide/write_docs.md @@ -0,0 +1,33 @@ +# Write related docs + +The most important thing after you implement your code in DMFF is to write the docs. Good documentation can help users use your module correctly and allow other maintainers to improve functions better. + +The documentation of DMFF use [MKDocs](https://www.mkdocs.org/) framework. + +## install MKDocs + +Before we start to write our docs, run the following command in termianl to install MKDocs: + +``` +pip install mkdocs +``` + +To support latex rendering, we also need markdown extension: + +``` +pip install pymdown-extensions +``` + +## Write your docs + +According to the existing document architecture, create a new markdown file in the appropriate directory. Write it! If you need to insert picture, upload the picture in `assets` directory, and use `![_placeholder](relative/path/to/this/file)` syntax to insert picture. + +## Preview you docs + +MkDocs comes with a built-in dev-server that lets you preview your documentation as you work on it. Make sure you're in the same directory as the `mkdocs.yml` configuration file, and then start the server by running the mkdocs serve command: + +``` +mkdocs serve +``` + +Open up http://127.0.0.1:8000/ in your browser, and you'll see the default home page being displayed. The dev-server also supports auto-reloading, and will rebuild your documentation whenever anything in the configuration file, documentation directory, or theme directory changes. \ No newline at end of file diff --git a/docs/index.md b/docs/index.md index 000ea3455..5406747fb 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,17 +1,27 @@ -# Welcome to MkDocs +# DMFF -For full documentation visit [mkdocs.org](https://www.mkdocs.org). +**DMFF** (**D**ifferentiable **M**olecular **F**orce **F**ield) is a Jax-based python package that provides a full differentiable implementation of molecular force field models. This project aims to establish an extensible codebase to minimize the efforts in force field parameterization, and to ease the force and virial tensor evaluations for advanced complicated potentials (e.g., polarizable models with geometry-dependent atomic parameters). Currently, this project mainly focuses on the molecular systems such as: water, biological macromolecules (peptides, proteins, nucleic acids), organic polymers, and small organic molecules (organic electrolyte, drug-like molecules) etc. We support both the conventional point charge models (OPLS and AMBER like) and multipolar polarizable models (AMOEBA and MPID like). The entire project is backed by the XLA technique in JAX, thus can be "jitted" and run in GPU devices much more efficiently compared to normal python codes. -## Commands +The behavior of organic molecular systems (e.g., protein folding, polymer structure, etc.) is often determined by a complex effect of many different types of interactions. The existing organic molecular force fields are mainly empirically fitted and their performance relies heavily on error cancellation. Therefore, the transferability and the prediction power of these force fields are insufficient. For new molecules, the parameter fitting process requires essential manual intervention and can be quite cumbersome. In order to automate the parametrization process and increase the robustness of the model, it is necessary to apply modern AI techniques in conventional force field development. This project serves for this purpose by utilizing the automatic differentiable programming technique to develop a codebase, which allows a more convenient incorporation of modern AI optimization techniques. It also helps the realization of many exciting functions including (but not limited to): hybrid machine learning/force field models and parameter optimization based on trajectory. -* `mkdocs new [dir-name]` - Create a new project. -* `mkdocs serve` - Start the live-reloading docs server. -* `mkdocs build` - Build the documentation site. -* `mkdocs -h` - Print help message and exit. +## User Guide -## Project layout ++ [1. Introduction](user_guide/introduction.md) ++ [2. Installation](user_guide/installation.md) ++ [3. Compute energy and forces](user_guide/compute.md) ++ [4. Compute gradients with auto differentiable framework](user_guide/auto_diff.md) ++ [5. Theories](user_guide/theory.md) ++ [6. Introduction to force field xml files](user_guide/xml_spec.md) - mkdocs.yml # The configuration file. - docs/ - index.md # The documentation homepage. - ... # Other markdown pages, images and other files. +## Developer Guide ++ [1. Introduction](dev_guide/introduction.md) ++ [2. Architecture](dev_guide/arch.md) ++ [3. Convention](dev_guide/convention.md) + +## Modules ++ [1. ADMP](modules/admp.md) + + +## Support and Contribution + +Please visit our repository on [GitHub](https://github.com/deepmodeling/DMFF) for the library source code. Any issues or bugs may be reported at our issue tracker. All contributions to DMFF are welcomed via pull requests! diff --git a/docs/modules/admp.md b/docs/modules/admp.md new file mode 100644 index 000000000..02cab44d7 --- /dev/null +++ b/docs/modules/admp.md @@ -0,0 +1,5 @@ +# ADMP + +Automatic Differentiable Multipolar Polarizable (ADMP) force field calculator. + +[Introduction](../admp/readme.md) diff --git a/docs/user_guide/installation.md b/docs/user_guide/installation.md new file mode 100644 index 000000000..4798144e7 --- /dev/null +++ b/docs/user_guide/installation.md @@ -0,0 +1,41 @@ +# 2. Installation +## 2.1 Install Dependencies ++ Install [jax](https://github.com/google/jax) (select the correct cuda version, see more details in the Jax installation guide): +```bash +pip install jax[cuda11_cudnn82] -f https://storage.googleapis.com/jax-releases/jax_releases.html +``` ++ Install [jax-md](https://github.com/google/jax-md): +```bash +pip install jax-md +``` ++ Install [OpenMM](https://openmm.org/): +```bash +conda install -c conda-forge openmm +``` +## 2.2 Install DMFF from Source Code +One can download the DMFF source code from github: +```bash +git clone https://github.com/deepmodeling/DMFF.git +``` +Then you may install DMFF using `pip`: +```bash +cd dmff +pip install . --user +``` + +## 2.3 Test Installation +To test if DMFF is correctly installed, you can run the following commands in an interactive python session: +```python +>>> import dmff +>>> import dmff.admp +``` + +You can also run the example scripts to test whether DMFF is installed correctly. +```bash +cd ./examples/water_1024 +python ./run_admp.py + +cd ./examples/water_pol_1024 +python ./run_admp.py +``` +Note that the scripts will run slower than expect if `DO_JIT = True` in `dmff/settings.py`. This is because the programm will do the jit compilation when a function is invoked in the first time. diff --git a/docs/user_guide/introduction.md b/docs/user_guide/introduction.md new file mode 100644 index 000000000..fc8033ec1 --- /dev/null +++ b/docs/user_guide/introduction.md @@ -0,0 +1,17 @@ +# 1. Introduction + +In this user guide, you will learn how to: + +- Install DMFF +- Compute energy and force +- Auto differentiate potential functions +- Couple DMFF with MD engine + +The first thing you should know is that DMFF is not an actual force field (such as OPLS or AMBER), but a differentiable implementation of various force field (or "potential") functional forms. It contains many modules: + +- ADMP module: Automatic Differentiable Multipolar Polarizable potential (MPID like potentials) +- Classical module: implements classical force fields (OPLS or GAFF like potentials) +- SGNN module: implements subgragh neural network model for intramolecular interactions + +Each module implements a particular form of force field, which takes a unique set of input parameters, usually provided in a XML file. With DMFF, one can easily compute the energy as well as energy gradients including: forces, virial tensors, and gradients with respect to force field parameters etc. + diff --git a/docs/user_guide/multipole_pme.md b/docs/user_guide/multipole_pme.md new file mode 100644 index 000000000..d207f2914 --- /dev/null +++ b/docs/user_guide/multipole_pme.md @@ -0,0 +1,605 @@ +# Notes on Multipolar Ewald + + + +## References: + +* [Anthony's book](https://oxford.universitypressscholarship.com/view/10.1093/acprof:oso/9780199672394.001.0001/acprof-9780199672394) +* The Multipolar Ewald paper in JCTC: ***J. Chem. Theory Comput.*** 2015, **11**, 2, 436–450 [Link](https://pubs.acs.org/doi/abs/10.1021/ct5007983) +* The `convert_mom_to_xml.py` code, which converts the CAMCASP moments (global harmonics) to OpenMM moments (local cartesian). +* The dispersion Ewald/PME: [Link](https://aip.scitation.org/doi/pdf/10.1063/1.470117) + + + +## Relevant mathematics + +### Cartesian Tensors + +Multipole operators in Cartesian Tensors: + +$$ +\displaylines{ +\begin{align} +&\hat{q} = 1 \\ +&\hat{\mu}_{\alpha} = r_{\alpha} \\ +&\hat{\Theta}_{\alpha\beta}=\frac{3}{2}\left(r_{\alpha}r_{\beta}-\frac{1}{3}r^2\delta_{\alpha\beta}\right) +\end{align} +} +$$ + +Interaction Hamiltonian (up to quadrupole): + +$$ +\displaylines{ +\hat{H}' = T q^{A} q^{B}+T_{\alpha}\left(q^{A} \hat{\mu}_{\alpha}^{B}-\hat{\mu}_{\alpha}^{A} q^{B}\right)+T_{\alpha \beta}\left(\frac{1}{3} q^{A} \hat{\Theta}_{\alpha \beta}^{B}-\hat{\mu}_{\alpha}^{A} \hat{\mu}_{\beta}^{B}+\frac{1}{3} \hat{\Theta}_{\alpha \beta}^{A} q^{B}\right) \\ +-\frac{1}{3} T_{\alpha \beta \gamma}\left(\hat{\mu}_{\alpha}^{A} \hat{\Theta}_{\beta \gamma}^{B}-\hat{\Theta}_{\alpha \beta}^{A} \hat{\mu}_{\gamma}^{B}\right) +\frac{1}{9}T_{\alpha \beta \gamma \delta}\hat{\Theta}_{\alpha \beta}^{A} \hat{\Theta}_{\gamma \delta}^{B} + \cdots +} +$$ + +And the corresponding interaction tensors are: + +$$ +T_{\alpha \beta \ldots v}^{(n)}=\frac{1}{4 \pi \epsilon_{0}} \nabla_{\alpha} \nabla_{\beta} \ldots \nabla_{v} \frac{1}{R} +$$ + +The first few terms are listed in the page 4 of Anthony's book (up to quad-quad interactions). + +### Spherical Tensors + +The often seen spherical harmonics include: + +* The original spherical harmonics (**adopting the convention without $(-1)^m$**): + +$$ +Y_{lm}(\theta, \phi) = \sqrt{\frac{2 l+1}{4 \pi} \frac{(l-m) !}{(l+m) !}} P_{l}^{m}(\cos \theta) e^{i m \phi} +$$ + + with $P_l^m$ being the associated Legendre Polynomials: +$$ +P_{\ell}^{m}(x)=\frac{(-1)^{m}}{2^{\ell} \ell !}\left(1-x^{2}\right)^{m / 2} \frac{d^{\ell+m}}{d x^{\ell+m}}\left(x^{2}-1\right)^{\ell} +$$ + Note we have: +$$ +P_{\ell}^{-m}(x)=(-1)^{m} \frac{(\ell-m) !}{(\ell+m) !} P_{\ell}^{m}(x) +$$ + The exact form of these polynomials can be found in: + + [Associated_Legendre_polynomials](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials) + + The lowest few orders of them are: + +$$ +\displaylines{ +\begin{aligned} +&P_{0}^{0}(x)=1 \\ +&P_{1}^{0}(x)=x \\ +&P_{1}^{1}(x)=-\left(1-x^{2}\right)^{1 / 2} \\ +&P_{2}^{0}(x)=\frac{1}{2}\left(3 x^{2}-1\right) \\ +&P_{2}^{1}(x)=-3 x\left(1-x^{2}\right)^{1 / 2} \\ +&P_{2}^{2}(x)=3\left(1-x^{2}\right) +\end{aligned} +} +$$ + +* We also use the renormalized spherical harmonics: + +$$ +\displaylines{ +C_{l m}(\theta, \varphi)=\sqrt{\frac{4 \pi}{2 l+1}} Y_{l m}(\theta, \varphi) \\ + = \sqrt{\frac{(l-m) !}{(l+m) !}} P_{l}^{m}(\cos \theta) e^{i m \phi} \\ + =\epsilon_m\sqrt{\frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|}(\cos\theta)e^{im\phi} +} +$$ + +**Note that I believe the equation B.1.3 in Pg. 272 of Anthony's book is not quite right** (missing an absolute value in $P_l^m$). Definition of $\epsilon_m$ can be found in Pg. 272. + +* Now we can define the regular and irregular spherical harmonics: + +$$ +\displaylines{ +\begin{aligned} +R_{l m}(\vec{r}) &=r^{l} C_{l m}(\theta, \varphi) \\ +I_{l m}(\vec{r}) &=r^{-l-1} C_{l m}(\theta, \varphi) +\end{aligned} +} +$$ + +In particular: + +$$ +R_{lm}(\vec{r}) = r^l \cdot \sqrt{\frac{(l-m) !}{(l+m) !}} P_{l}^{m}(\cos \theta) e^{i m \phi} +$$ + +From this we can define the real-valued regular spherical harmonics (for $m>0$): + +$$ +\displaylines{ +\begin{aligned} +R_{l m c} &=\sqrt{\frac{1}{2}}\left[(-1)^{m} R_{l m}+R_{l,-m}\right] \\ +\mathrm{i} R_{l m s} &=\sqrt{\frac{1}{2}}\left[(-1)^{m} R_{l m}-R_{l,-m}\right] +\end{aligned} +} +$$ + +According to my derivation, that is: + +$$ +\displaylines{ +\begin{cases} +R_{l0} = r^l\cdot C_{l0} = r^l\cdot P_l^{0}(\cos\theta)\\ +R_{lmc} = (-1)^m \cdot \sqrt{2} \cdot r^l \sqrt{\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)\cos(m\phi) \\ +R_{lms} = (-1)^m \cdot \sqrt{2} \cdot r^l \sqrt{\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)\sin(m\phi) +\end{cases} +} +$$ + +In Anthony's book, they usually use Greek letters ($\mu,\nu,\kappa$ etc.​) or letter $t,u$ to represent $1c, 1s, 2c, 2s$​ etc. These are the projector functions to compute all the moments. + +In the real-valued regular spherical harmonics representation, the interaction energy is: + +$$ +\hat{H}' = \sum_{at,bu} {\hat{Q}_t^a T_{tu}^{ab}\hat{Q}_u^b} +$$ + +The interaction tensor $T_{tu}^{AB}$ can be found in Appendix F of Anthony's book (pg. 291). + +* A couple of things regarding paper: ***J. Chem. Theory Comput.*** 2015, **11**, 2, 436–450 (https://pubs.acs.org/doi/abs/10.1021/ct5007983) + +In this paper, the definition (***in their notation***) of: + +$$ +\displaylines{ +C_{l\mu} = (-1)^{\mu}\sqrt{\left(2-\delta_{\mu, 0}\right)(l+\mu) !(l-\mu) !}\cdot R_{l\mu} \\ += +\begin{cases} +(-1)^{m}\sqrt{\left(2-\delta_{m, 0}\right)(l+m) !(l-m) !} \cdot R_{l mc}(\vec{r}) & \mu \geq 0 \\ +(-1)^{-m}\sqrt{2(l-m) !(l+m) !} \cdot R_{lms}(\vec{r}) & \mu<0 +\end{cases} \\ +(\text{Here, we assumes $m=|\mu|$}) \\ += +\begin{cases} +r^l\cdot P_l^0(\cos\theta) & \mu=0 \\ +(-1)^m\cdot\sqrt{2}\cdot r^l\sqrt{\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)\cos(m\phi) & \mu\gt 0 \\ +(-1)^m\cdot\sqrt{2}\cdot r^l\sqrt{\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)\sin(m\phi) & \mu\lt 0 +\end{cases} +} +$$ + +This is identical to Eqn. 14. **Therefore, the $C_{l\mu}$ in the JCTC paper is identical to the $R_\mu$ in Anthony's book**. + +And the interaction function also checks out! So we have a compact form for the interaction tensor in table F.1 of the book + +$$ +T_{tu}^{ab} = \frac{R_t(\nabla_a)}{(2l_t-1)!!}\frac{R_u(\nabla_b)}{(2l_u-1)!!}\frac{1}{R_{ab}} +$$ + +This is verified manually by computing the $T_{21c,11c}=T_{xz, x} = \sqrt{3}\frac{1}{R^4}\left(\hat{z}-5\hat{x}^2\hat{z}\right)$​ tensor. + + + +## Rotations + +Assume the local frame matrix for a particular atom is: + +$$ +\displaylines{ +\mathbf{R} = +\left[ +\matrix{ +\mathbf{r_1^T} \\ +\mathbf{r_2^T} \\ +\mathbf{r_3^T} \\ +}\right] = \left[ +\matrix{ +x_1 & y_1 & z_1 \\ +x_2 & y_2 & z_2 \\ +x_3 & y_3 & z_3 +}\right] +} +$$ + +Then for any vector, the local ($\vec{r}'$)-to-global (\vec{r}) transform is: + +$$ +\vec{r}=\mathbf{R}^T\cdot\vec{r}' +$$ + +Following the convention in the `convert_mom_to_xml.py` code, assuming the order for spherical harmonics is: $Q_{00}, Q_{10}, Q_{11c}, Q_{11s}, Q_{20}, Q_{21c}, Q_{21s}, Q_{22c}, Q_{22s}$​​ + +* **Dipole rotation** + +For dipole, the harmonics-to-cartesian conversion is: + +$$ +\mathbf{C}_1=\left[ +\matrix{ +0 & 1 & 0 \\ +0 & 0 & 1 \\ +1 & 0 & 0 +} +\right] +$$ + +So the **local-harmonics-to-global-harmonics** conversion for dipole moments is: + +$$ +\mathbf{R}_{1}^{l-g}=\mathbf{C}_1^T\mathbf{R}^T\mathbf{C}_1 +$$ + +* **Quadrupole rotation** + +First of all, the cartesian moments are: + +$$ +\displaylines{ +\mathbf{\Theta}= +\int {d\vec{r}\cdot\rho(\vec{r})\cdot\left(\frac{3}{2}\left[\matrix{ +xx & xy & xz \\ +yx & yy & yz \\ +zx & zy & zz \\ +}\right] - \frac{1}{2}r^2\mathbf{I} \right)} \\ += \int {d\vec{r}\cdot\rho(\vec{r})\cdot\left(\frac{3}{2}\vec{r}\vec{r}^T -\frac{1}{2}(\vec{r}^T\cdot\vec{r})\mathbf{I}\right)} +} +$$ + +Then obviously, the **local-to-global** conversion is: + +$$ +\mathbf{\Theta} = \mathbf{R}^T\mathbf{\Theta}'\mathbf{R} +$$ + +The conversion between $\vec{\Theta}=[xx,yy,zz,xy,xz,yz]$​​​​ and harmonic moments ($\vec{Q}=[Q_{20}, Q_{21c}, Q_{21s}, Q_{22c}, Q_{22s}]$​​​​​) is: + +$$ +\displaylines{ +\mathbf{C}_2^{h-c} = \left[\matrix{ +-\frac{1}{2} & 0 & 0 & \frac{\sqrt{3}}{2} & 0 \\ +-\frac{1}{2} & 0 & 0 & -\frac{\sqrt{3}}{2} & 0 \\ +1 & 0 & 0 & 0 & 0 \\ +0 & 0 & 0 & 0 & \frac{\sqrt{3}}{2} \\ +0 & \frac{\sqrt{3}}{2} & 0 & 0 & 0 \\ +0 & 0 & \frac{\sqrt{3}}{2} & 0 & 0 +}\right] +} +$$ + +$$ +\vec{\Theta} = \mathbf{C}_2^{h-c} \vec{Q} +$$ + +And the inversed conversion is: + +$$ +\displaylines{ +\mathbf{C}_2^{c-h} = \left[\matrix{ +0 & 0 & 1 & 0 & 0 & 0 \\ +0 & 0 & 0 & 0 & \frac{2}{\sqrt{3}} & 0 \\ +0 & 0 & 0 & 0 & 0 & \frac{2}{\sqrt{3}} \\ +\frac{1}{\sqrt{3}} & -\frac{1}{\sqrt{3}} & 0 & 0 & 0 & 0 \\ +0 & 0 & 0 & \frac{2}{\sqrt{3}} & 0 & 0 \\ +}\right] +} +$$ + +$$ +\vec{Q} = \mathbf{C}_2^{c-h}\vec{\Theta} +$$ + +So the rotation procedure from **local-to-global** rotation for (spherical harmonics) quadrupole moments is: + +$$ +\displaylines{ +\vec{\Theta}' = \mathbf{C}_2^{h-c} \vec{Q}' \\ +\mathbf{\Theta} = \mathbf{R}^T\mathbf{\Theta}'\mathbf{R} \\ +\vec{Q} = \mathbf{C}_2^{c-h}\vec{\Theta} +} +$$ + +And vice-versa for the invert process. + +Using `sympy`, we can derive the actual linear transformation matrix, which is coded in `test_rot_quad.py`. The derivation is done using `derive_rot_mat.py` + + + + + +## Multipolar Ewald + +* **The reciprocal space part:** + + The structure factor is: + + $$ + S(\vec{k}) = \sum_{a,t} Q_t^a\frac{R_t(\vec{k})(-i)^{l_t}}{(2l-1)!!}e^{-i\vec{k}\cdot\vec{R}_a} + $$ + + The reciprocal space energy is: + + $$ + E_{recip} = \sum_{k\neq 0} {\frac{2\pi}{V k^2}\exp\left(-\frac{k^2}{4\kappa^2}\right)\left|S(\vec{k})\right|^2} + $$ + +* The real space part: + +$$ +\displaylines{ +E_{real} = \sum_{\substack{a + + + + + + + + + + + + + + + + + + +``` + +Where "-C" indicates an external bond between the "C" atom and the **previous** residue. During typification, the code will try to match all atoms in the topology template with the actual atoms in the real structure. Once a match is successful, the matched atom will be bonded accordingly. If the match fails, this template atom will be skipped. Therefore, the actual number of bonds in the matched structure can be less than the number of bonds defined in the template. + +The XML file registration method is as follows: + +``` python +try: + import openmm.app as app +except: + import simtk.openmm.app as app + +app.Topology.loadBondDefinations("residues.xml") # register residue topology + +# Create topology and add atoms and residues to it, which is automatically performed when reading PDB +top = app.Topology() +... +top.createStandardBonds() # Connect keys according to template files +``` + +After this process, the bonding topologies are constructed in the matched residue, but the force field parameters are not yet assigned. It should be noted that disulfide bonds are not registered in this step. The OpenMM topology class will look for SG atoms in Cys that are not connected to Hg, and connect SG atom pairs within 0.3nm as disulfide bonds. + +## Force field Parameter File + +After all bonds are constructed using the topology file, the force field parameters will be assigned using the force field parameter file. + +The force field parameter file is as follows: +``` xml + + + + + + + + + + + + + + + + + + + + + + + + + + + +``` +This file can be further divided into the residue part and the force field part. + +### Residue Part +``` xml + + + + + + + + + + + + ... + +``` +The `` node of the residue part defines all the atoms involved in the residue and some paramemters per atom, which can be called by the force field part on demand. The `` node defines the bonding information of the residue. The information contained in this part is different from that in the topology file discussed above. Take ALA as an example, we usually have at least three states for ALA, N-end, C-end and in-chain. The corresponding parameter templates in the force field file are as follows: + +``` xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + <\displaylines{Bond atomName1="CA" atomName2="HA"/> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +``` + +In this example, the atom numbers and the bonding configurations of ALA, CALA and NALA are different. When matching each ALA, OpenMM will try to match CALA, NALA, and ALA separately. It will compare each parameter template with the topology of the residue, and select the one with the right number of atoms, element composition, and bonding configurations as the matched template. The parameter template contains atom type and class information, which are then used to assign force field parameters. + + +### Forcefield Part +```xml + + + ... + + + + + + + + + + + + + + + + +``` + +The `` node defines all atom types. The `type` label of each atom in the residue part will match the `name` label of each child node of ``. For each atom type, it also defines a `class` tag for different matching scenarios. The `name` tag of different `` must be different, but the `class` tag can be the same. + +The `<*Force>` node defines the matching rule of a potential function. For example, `` defines harmonic bond, and the `` node defines intermolecular interaction. For more information about each force, the readers are referred to this document: [details](http://docs.openmm.org/latest/userguide/application/05_creating_ffs.html#writing-the-xml-file). + +In the matching process, OpenMM will iterate all atoms, bonds, angles, dihedrals and impropers, and add all matched entries to the total potential function. Matching can be carried out according to the `type` tag, corresponding to the `name` of each atom defined in ``; It can also be based on the `class` tag, corresponding to the `class` of each atom in ``. This design is applicable to the situation when there are many types of atoms but they are roughly the same. For example, there are few kinds of LJ parameters in small molecular force field, but there are many kinds of intramolecular force parameters. We can even create a separate type for a specific small molecule to define the intra molecular interaction, but it belongs to the same class on LJ, so as to achieve the effect that the small molecule parameters can be tuned and do not affect each other. diff --git a/mkdocs.yml b/mkdocs.yml index 14a7a7497..6833919a7 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -1,9 +1,32 @@ site_name: DMFF nav: - - Home: 'index.md' - - About: 'about.md' - - 'User Guide': - - 'tutorial': 'tutorial.md' - - 'Modules': - - 'ADMP': + - Home: index.md + - User Guide: + - 1. Introduction: user_guide/introduction.md + - 2. Installation: user_guide/installation.md + - 3. Compute energy and force: user_guide/compute.md + - 4. Auto-diff: user_guide/auto_diff.md + - 5. Theory: user_guide/theory.md + - 6. couple with MD: user_guide/couple.md + - 7. Introduction to force field xml file: user_guide/xml_spec.md + + - Developer Guide: + - 1. Introduction: dev_guide/introduction.md + - 2. Architecture: dev_guide/arch.md + - 3. Convention: dev_guide/convention.md + - Modules: + - ADMP: + - Introduction: admp/readme.md + - Frontends: admp/frontend.md + - About: about.md + theme: readthedocs + +markdown_extensions: + - pymdownx.arithmatex: + generic: true + +extra_javascript: + - javascripts/mathjax.js + - https://polyfill.io/v3/polyfill.min.js?features=es6 + - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js