Merge pull request #135 from Cascella-Group-UiO/cbt

Electrostatic and dihedrals

.gitignore

@@ -14,3 +14,8 @@ coverage.xml
 __pycache__
 RUN/
 *.so
+test-xinmeng-electrostatic/*
+*.h5
+!examples/*.h5
+!test-xinmeng-electrostatic/run-dppc/
+#!test-xinmeng-electrostatic/run-sds/

README.md

@@ -1,18 +1,37 @@
 HyMD testing and development · [![License: GPL v3](https://img.shields.io/badge/License-LGPLv3-blue.svg)](https://www.gnu.org/licenses/lgpl-3.0.html) ![build](https://github.com/Cascella-Group-UiO/HyMD-2021/workflows/build/badge.svg)
 ---------
-Compile FORTRAN modules:
+HyMD is a software that can run coarse-grained molecular dynamics simulations
+using the hybrid-particle field model approach, introduced initially in [1].
+In HyMD we implement the formulation presented in [2].
+
+## Setup
+First compile the FORTRAN modules:
 ```bash
 > cd hymd/
 > make clean
 > make
 > cd ..
 ```
 
-Run a simple example simulation with
+Check out all the available options with
 ```bash
-> mpirun -n 4 python3 hymd/main.py config.toml dppc.h5 --verbose -logfile log.txt
+> python3 hymd/main.py --help
 ```
 
-Notes:
-- Changed numpy data types (Float32, Int32 -> float32, int32) to run in local computer
-- TODO: Fix for running monoatomic (or non molecular) systems
+Start a simple example simulation with:
+```bash
+> mpirun -n 4 python3 hymd/main.py peptide.toml peptide.h5 -v
+```
+
+## References
+[1] Milano, G. & Kawakatsu, T. Hybrid particle-field molecular dynamics
+simulations for dense polymer systems Journal of Chemical Physics, American
+Institute of Physics, 2009, 130, 214106 
+
+[2] Bore, S. L. & Cascella, M. Hamiltonian and alias-free hybrid
+particle--field molecular dynamics The Journal of Chemical Physics, AIP
+Publishing LLC, 2020, 153, 094106 
+
+
+
+

config.toml

@@ -1,57 +1,82 @@
-[meta]
+# [meta]
 # Name of the simulation. May be ommitted.
-name = "DPPC bilayer with ML interaction parameters"
+# name = "PRWG lipopeptide micelle in water"
 # Tags classifying the simulation. May be ommitted.
-tags = ["bilayer", "solvent", "DPPC"]
+# tags = ["lipopeptides", "PRWG"]
+
+# Here's the topology of a PRWG molecule
+#
+#                   WSC2 -- WSC3
+#                   |       |
+#    PSC(+)        WSC1 -- WSC4
+#    |             |
+#    PBB -- RBB -- WBB -- GBB -- C -- C -- C -- C 
+#           |
+#           RSC -- RSC(+)
+#
 
 [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
+# ommited here.
+n_particles = 11876
+# specify number of ghost particles, used for dipole point charges in the CBT dihedral potential (default = 0)
+# n_ghots = 0
 # 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
+# max_molecule_size = 15
 
 [simulation]
-# Number of total time steps in the simulation in [picoseconds].
-n_steps = 25000
+# Number of total time steps in the simulation.
+n_steps = 200
 # Frequency of trajectory/energy file output in time steps.
-n_print = 5000
+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
+n_flush = 10
 # Time step used in the simulation in [picoseconds].
-time_step = 0.3
+time_step = 0.03
 # Simulation box size in [nanometers].
-box_size = [13.0, 13.0, 14.0]
+box_size = [10.73911, 10.76184, 11.30546]
+# box_size = [5.0, 5.0, 5.0] # for single PRWG molecule test
 # 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
+respa_inner = 5
 # 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
+# domain_decomposition = false
+# Remove linear center of mass momentum from the system before integration
+# starts.
+cancel_com_momentum = true
 # Starting temperature to generate before simulation begins in [kelvin]. Ommit
 # or set to 'false' to not change the temperature before starting.
-start_temperature = 323
+start_temperature = 300
 # 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
+target_temperature = 300
 # Thermostat collision frequency in [1/picoseconds].
-tau = 0.1
+tau = 1
+# Couple groups of particles species to individual different thermostats.
+thermostat_coupling_groups = [
+  ["PBB", "PSC", "RBB", "RSC", "WBB", "WSC1", "WSC2", "WSC3", "WSC4", "GBB", "C"], 
+  ["W", "CL"]
+]
 # The energy functional W[phi] to use. Options:
 #    "SquaredPhi":      φ² / 2κφ₀,
 #    "DefaultNoChi":   (φ - φ₀)² / 2κφ₀
 #    "DefaultWithChi": (φ - φ₀)² / 2κφ₀ + Σ χφφ' / φ₀
 # Subclass Hamiltonian to create a new energy functional.
 hamiltonian = 'DefaultWithChi'
+coulombtype = 'PIC_Spectral'
+dielectric_const = 80.0
 
 [field]
 # Particle-mesh grid size, either a single integer or an array of 3 integers
@@ -64,39 +89,9 @@ mesh_size = [24, 24, 24]
 # [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
+# This value is a characteristic 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],
-]
+chi = []

examples/dppc.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 = 25000
+# Frequency of trajectory/energy file output in time steps.
+n_print = 5000
+# 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]],
+]