Merge pull request #206 from samiransen23/main

Release 2.0.0

.gitignore

@@ -1,7 +1,5 @@
 build
 _build
-_static
-_templates
 hymd.egg-info
 *.hypothesis
 *.mypy_cache
@@ -17,6 +15,6 @@ RUN/
 *.so
 test-xinmeng-electrostatic/*
 *.h5
-!examples/*.h5
+examples/
 !test-xinmeng-electrostatic/run-dppc/
 #!test-xinmeng-electrostatic/run-sds/

README.md

@@ -7,9 +7,12 @@
 ---------
 **HylleraasMD** (HyMD) is a massively parallel Python package for hybrid particle-field molecular dynamics (hPF-MD) simulations of coarse-grained bio- and soft-matter systems.
 
-HyMD can run canonical hPF-MD simulations [[1]](#1), or filtered density Hamiltonian hPF (HhPF-MD) simulations [[2]](#2), with or without explicit PME electrostatic interactions [[3]](#3). It includes all standard intramolecular interactions, including stretching, bending, torsional, and combined bending-dihedral potentials. Additionally, topological reconstruction of permanent peptide chain backbone dipoles is possible for accurate recreation of protein conformational dynamics [[4]](#4). Martini style elastic networks (ElNeDyn) [[5]](#5) are also supported.
+HyMD can run canonical hPF-MD simulations, or filtered density Hamiltonian hPF (HhPF-MD) simulations [[1]](#1),[[2]](#2),[[3]](#3) with or without explicit PME electrostatic interactions. It includes all standard intramolecular interactions,
+including stretching, bending, torsional, and combined bending-dihedral potentials. Additionally, topological reconstruction of permanent peptide chain backbone dipoles is possible for accurate recreation of protein conformational dynamics.
+It can run simulations in constant energy (NVE), constant volume (NVT) [[1]](#1) or constant pressure (NPT) conditions [[4]](#4).
 
-HyMD uses the [pmesh](github.com/rainwoodman/pmesh) library for particle-mesh operations, with the PPFT [[6]](#6) backend for FFTs through the [pfft-python bindings](github.com/rainwoodman/pfft-python). File IO is done via HDF5 formats to allow MPI parallel reads.
+HyMD uses the [pmesh](github.com/rainwoodman/pmesh) library for particle-mesh operations, with the PPFT [[5]](#5) backend for FFTs through the [pfft-python bindings](github.com/rainwoodman/pfft-python).
+File IO is done via HDF5 formats to allow MPI parallel reads.
 
 ## User Guide
 Detailed installation and user guide, together with comprehensive example simulations are located in the [HylleraasMD documentation](https://cascella-group-uio.github.io/HyMD/index.html).
@@ -73,23 +76,24 @@ chmod +x pytest-mpi
 pytest-mpi -oo -n 2 -ns
 ```
 
+## Please cite our work
+If you use HyMD for your purposes, please cite the appropriate references from the section below.  
+If you cannot cite all, the fundamental works to be cited are [[1]](#1) and [[4]](#4). 
+
 ---------
 
 ### References
 <a id="1">[1]</a>
-Milano, G.; Kawakatsu, T. Hybrid particle-field molecular dynamics simulations for dense polymer systems. J. Chem. Phys. **2009**, 130, 214106.
+Ledum, M.; Sen, S.; Li, X.; Carrer, M.; Feng Y.; Cascella, M.; Bore, S. L. HylleraasMD: A Domain Decomposition-Based Hybrid Particle-Field Software for Multi-Scale Simulations of Soft Matter. ChemRxiv 2021
 
 <a id="2">[2]</a>
-Bore, S. L.; Cascella, M. Hamiltonian and alias-free hybrid particle–field molecular dynamics. J. Chem. Phys. **2020**, 153, 094106.
+Ledum, M.; Carrer, M.; Sen, S.; Li, X.; Cascella, M.; Bore, S. L. HyMD: Massively parallel hybrid particle-field molecular dynamics in Python.
 
 <a id="3">[3]</a>
-Kolli, H. B.; De Nicola, A.; Bore, S. L.; Schäfer, K.; Diezemann, G.; Gauss, J.; Kawakatsu, T.;Lu, Z.-Y.; Zhu, Y.-L.; Milano, G.; Cascella, M. Hybrid Particle-Field Molecular DynamicsSimulations of Charged Amphiphiles in an Aqueous Environment. J. Chem. Theory Comput. **2018**, 14, 4928–4937.
+Bore, S. L.; Cascella, M. Hamiltonian and alias-free hybrid particle–field molecular dynamics. J. Chem. Phys. 2020, 153, 094106.
 
 <a id="4">[4]</a>
-Bore, S. L.; Milano, G.; Cascella, M. Hybrid Particle-Field Model for Conformational Dynamics of Peptide Chains. J. Chem. Theory Comput. **2018**, 14, 1120–1130.
+Sen, S.; Ledum, M.; Bore, S. L.; Cascella, M. Soft Matter under Pressure: Pushing Particle–Field Molecular Dynamics to the Isobaric Ensemble. ChemRxiv 2023
 
 <a id="5">[5]</a>
-Periole, X.; Cavalli, M.; Marrink, S. J.; Ceruso, M. A. Combining an elastic network with a coarse-grained molecular force field: structure, dynamics, and intermolecular recognition. J. Chem. Theory Comput. **2009**, 5.9, 2531-2543.
-
-<a id="6">[6]</a>
-Pippig, M. PFFT: An extension of FFTW to massively parallel architectures. SIAM J. Sci. Comput. **2013**, 35, C213–C236.
+Pippig, M. PFFT: An extension of FFTW to massively parallel architectures. SIAM J. Sci. Comput. 2013, 35, C213–C236.

config.toml

@@ -0,0 +1,102 @@
+[meta]
+# Name of the simulation. May be ommitted.
+name = "DPPC bilayer with ML interaction parameters"
+# Tags classifying the simulation. May be ommitted.
+tags = ["bilayer", "solvent", "DPPC"]
+
+[particles]
+# Number of total particles in the simulation. If an input .hdf5 format file is
+# specified, the number of particles will be inferred from this and *may* be
+# ommited.
+n_particles = 20336
+# Mass of the particles in [g/mol]. All masses are assumed equal.
+mass = 72.0
+# Maximum number of particles per molecules present in the system. A default of
+# 200 is assumed, and this keyword may be ommitted for any system with smaller
+# molecules.
+max_molecule_size = 15
+
+[simulation]
+# Number of total time steps in the simulation in [picoseconds].
+n_steps = 1
+# Frequency of trajectory/energy file output in time steps.
+n_print = 1
+# Frequency of requesting that the HDF5 library flush the file output buffers
+# to disk after in number of n_print timesteps.
+n_flush = 5000
+# Time step used in the simulation in [picoseconds].
+time_step = 0.3
+# Simulation box size in [nanometers].
+box_size = [13.0, 13.0, 14.0]
+# Time integrator used in the simulation. Either "velocity-verlet" or "respa".
+# If "respa", specify also the number of small rRESPA time steps per large
+# time_step with the 'respa_inner' keyword.
+integrator = "respa"
+respa_inner = 10
+# Perform MPI rank domain decomposition every x time steps to (hopefully)
+# reduce the amount of neccessary communication between ranks in the pmesh
+# procedures. Ommit or set to 'false' or '0' to not perform any domain
+# decomposition.
+domain_decomposition = 1000
+# Starting temperature to generate before simulation begins in [kelvin]. Ommit
+# or set to 'false' to not change the temperature before starting.
+start_temperature = 323
+# Target temperature used in the velocity rescale thermostat in [kelvin]. Ommit
+# or set to 'false' to use no thermostat, i.e. a constant energy simulation.
+target_temperature = 323
+# Thermostat collision frequency in [1/picoseconds].
+tau = 0.1
+# The energy functional W[phi] to use. Options:
+#    "SquaredPhi":      φ² / 2κφ₀,
+#    "DefaultNoChi":   (φ - φ₀)² / 2κφ₀
+#    "DefaultWithChi": (φ - φ₀)² / 2κφ₀ + Σ χφφ' / φ₀
+# Subclass Hamiltonian to create a new energy functional.
+hamiltonian = 'DefaultWithChi'
+
+[field]
+# Particle-mesh grid size, either a single integer or an array of 3 integers
+# (number of grid points in each dimension). In order to guarantee consistency
+# and speed in the PFFT routines, the actual mesh grid will be changed to ensure
+# that each dimension of the 2d PFFT process grid divides each dimension of the
+# mesh grid size.
+mesh_size = [24, 24, 24]
+# Compressibility used in the relaxed incompressibility term of W(phi) in
+# [mol/kJ].
+kappa = 0.05
+# Standard deviation in the Gaussian filter (window function) in [nanometers].
+# This value is a characzteristic length scale for the size of the particles in
+# the simulation.
+sigma = 0.5
+# Interaction matrix, chi, ((atom name 1, atom name 2), (mixing energy in
+# [kJ/mol])).
+chi = [
+  [["C", "W"], [42.24]],
+  [["G", "C"], [10.47]],
+  [["N", "W"], [-3.77]],
+  [["G", "W"],  [4.53]],
+  [["N", "P"], [-9.34]],
+  [["P", "G"],  [8.04]],
+  [["N", "G"],  [1.97]],
+  [["P", "C"], [14.72]],
+  [["P", "W"], [-1.51]],
+  [["N", "C"], [13.56]],
+]
+
+[bonds]
+# Two-particle bonds, ((atom name 1, atom name 2), (equilibrium length in
+# [nanometers], bond strenght in [kJ/mol])). Note the two
+bonds = [
+  [["N", "P"], [0.47, 1250.0]],
+  [["P", "G"], [0.47, 1250.0]],
+  [["G", "G"], [0.37, 1250.0]],
+  [["G", "C"], [0.47, 1250.0]],
+  [["C", "C"], [0.47, 1250.0]],
+]
+# Three-particle angular bonds, ((atom name 1, atom name 2, atom name 3),
+# (equilibrium angle in [degrees], bond strenght in [kJ/mol])).
+angle_bonds = [
+  [["P", "G", "G"], [120.0, 25.0]],
+  [["P", "G", "C"], [180.0, 25.0]],
+  [["G", "C", "C"], [180.0, 25.0]],
+  [["C", "C", "C"], [180.0, 25.0]],
+]

docs/doc_pages/api.rst

@@ -19,6 +19,24 @@ Thermostat
    :undoc-members:
 
 
+Barostat
+--------
+
+Berendsen Barostat
+^^^^^^^^^^^^^^^^^^
+
+.. automodule:: hymd.barostat
+   :members: isotropic, semiisotropic
+   :undoc-members:
+
+SCR Barostat
+^^^^^^^^^^^^^^^^^^
+
+.. automodule:: hymd.barostat_bussi
+   :members: isotropic, semiisotropic
+   :undoc-members:
+
+
 Hamiltonian
 -----------
 
@@ -55,6 +73,14 @@ Force
    :undoc-members:
 
 
+Pressure
+--------
+
+.. automodule:: hymd.pressure
+   :members: comp_pressure
+   :undoc-members:
+
+
 Logger
 ------
 

docs/doc_pages/config_file.rst

@@ -192,3 +192,53 @@ Configuration keywords specifying electrostatics and electrostatic parameters.
    :code:`float` [**optional**, default: :code:`None`]
 
    Specifies the relative dielectric constant of the simulation medium which regulates the strength of the electrostatic interactions. When using helical propensity dihedrals, this keyword must be specified---even if electrostatics are not included with the :code:`coulombtype` keyword.
+
+Pressure keywords
+^^^^^^^^^^^^^^^^^
+Configuration keywords specifying pressure and barostat parameters.
+
+Simulation keywords
+=========================================
+
+:pressure:
+    :code:`boolean` [**optional**, default: :code:`False`]
+
+    Specifies whether or not to calculate total internal pressure vector.
+
+:barostat:
+    :code:`string` [**optional**, default: :code:`None`](options: :code:`isotropic` or :code:`semiisotropic`)
+
+    Specifies whether to apply pressure constraints equally in all 3 Cartesian directions (`isotropic`) or equally in xy and different in z (`semiisotropic`).
+
+:barostat_type:
+    :code:`string` [**optional**, default: :code:`berendsen`](options: :code:`berendsen` or :code:`scr`)
+
+    Specifies the type of barostat to use. :code:`berendsen` is more suitable for equilibration and :code:`scr` for equilibrium data collection.
+
+:n_b:
+    :code:`integer` [**optional**, default: :code:`1`]
+
+    Frequency of barostat call in number of outer rRESPA steps. 
+
+:tau_p:
+   :code:`float` [**optional**, default: :code:`10 tau` if :code:`tau < 0.1` else :code:`1.0`] {units: :math:`\text{ps}=10^{-12}~\text{s}`}
+
+   The time scale of the barostat coupling.
+
+:target_pressure:
+   :code:`array` [:code:`float`] [**optional**, default: :code:1] {units: `bar`}
+
+   Couples the system to an external pressure set to `target_pressure`.
+
+Field keywords
+=========================================
+
+:rho0:
+    :code:`float` [**optional**, default: :code:average density] {units: :math:`\text{nm}^{-3}`}
+
+    Intrinsic parameter corresponding to the specific volume of a coarse-grained particle. Typically: `8.66` 
+
+:a:
+    :code:`float` [:required:`required`] {units: :math:`\text{nm}^{-3}`}
+
+    Calibrated parameter to obtain the correct average density at the target temperature and pressure. Typically: `9.21`

docs/doc_pages/pressure.rst

@@ -0,0 +1,37 @@
+.. _pressure-label:
+
+Pressure
+########
+Internal pressure is calculated from internal energy according to
+
+.. math::
+
+    P_a = \frac{1}{\mathcal{V}} \left( 2T_a - \text{Vir}_a \right) \\
+    \text{Vir}_a = L_a \frac{\partial \mathcal{U}}{\partial L_a} \\
+    \mathcal{U} = \sum_{i=1}^M \mathcal{U}_0( \{ \mathbf{r}\}_i ) + W[\{ \tilde\phi \} ]
+
+where
+:math:`\mathcal{V}` is the simulation volume,
+:math:`{T_a}` is the kinetic energy
+and :math:`L_a` the length of the box in the Cartesian direction :math:`a`,
+Vir is the virial of the total interaction energy :math:`\mathcal{U}`.
+
+:math:`\mathcal{U}` comprises intramolecular bonded terms :math:`\mathcal{U}_0` (see :ref:`bonds-label` for details),
+and field terms :math:`W[\{ \tilde\phi \} ]` (see :ref:`theory-label` for details).
+
+Using the above expressions, the following form for internal pressure is obtained:
+
+.. math::
+
+    P_a = \frac{2 T_a}{\mathcal{V}} -\frac{L_a}{\mathcal{V}} \sum_{i=1}^N \frac{\partial \mathcal{U}_{0i}}{\partial L_a} + P^{(3)}_a \\
+
+.. P^{(3)}_a = - \frac{L_a}{\mathcal{V}} \frac{\partial W[\{ \tilde\phi \} ]}{\partial L_a} \\
+
+.. math::
+
+    P^{(3)}_a = \frac{1}{\mathcal{V}}\left ( -W[\{ \tilde\phi(\mathbf{r}) \}] + \int \sum_t \bar{V}_t(\mathbf{r})\tilde\phi_t(\mathbf{r})d\mathbf{r} 
+                + \int \sum_t \sigma^2\bar{V}_t(\mathbf{r})\nabla_a^2\tilde\phi_t(\mathbf{r}) d\mathbf{r} \right)
+
+where :math:`\bar{V}_t(\mathbf{r}) = \frac{\partial w(\{\tilde\phi\})}{\partial\tilde\phi_t}` 
+and :math:`σ` is a coarse-graining parameter (see :ref:`filtering-label` for details).
+Note that the above expression is obtained for a Gaussian filter which is the most natural choice in HhPF theory.

docs/doc_pages/theory.rst

@@ -107,6 +107,14 @@ The default form of the interaction energy functional in HyMD is
 
    W=\frac{1}{2\phi_0}\int\mathrm{d}\mathbf{r}\sum_{\text{i},\text{j}}\tilde\chi_{\text{i}-\text{j}}\tilde\phi_\text{i}(\mathbf{r})\tilde\phi_\text{j}(\mathbf{r}) + \frac{1}{2\kappa}\left(\sum_\text{k}\tilde\phi_\text{k}(\mathbf{r})-\phi_0\right)^2.
 
+In the case of constant pressure simulations (NPT), the interaction energy functional becomes
+
+.. math::
+
+   W=\frac{1}{2\rho_0}\int\mathrm{d}\mathbf{r}\sum_{\text{i},\text{j}}\tilde\chi_{\text{i}-\text{j}}\tilde\phi_\text{i}(\mathbf{r})\tilde\phi_\text{j}(\mathbf{r}) + \frac{1}{2\kappa}\left(\sum_\text{k}\tilde\phi_\text{k}(\mathbf{r})-a\right)^2.
+
+where, :math:`\rho_0= 1 / ν` is an intrinsic parameter corresponding to the specific volume :math:`(ν)` of a coarse-grained particle,  
+:math:`a` is a calibrated parameter to obtain the correct average density at the target temperature and pressure.
 See :ref:`functionals-label` for details.
 
 

docs/doc_pages/topology_input.rst

@@ -75,3 +75,10 @@ Optional datasets
    :code:`datatype` [:code:`float32` or :code:`float64`]
 
    The shape represents one particle charge per :code:`N` particles. The data type may be either four or eight byte reals.
+
+:/box:
+   :code:`array` shape [:code:`3`,]
+
+   :code:`datatype` [:code:`float32` or :code:`float64`]
+
+   The shape represents the initial box dimensions in Cartesian coordinates.

docs/index.rst

@@ -155,6 +155,7 @@ Indices and tables
    ./doc_pages/theory
    ./doc_pages/intramolecular_bonds
    ./doc_pages/electrostatics
+   ./doc_pages/pressure
 
 .. toctree::
    :maxdepth: 2