---
# System prepended metadata

title: HyMD - Generalized Possion Equation

---

# HyMD - Generalized Possion Equation 

### Status and in progress
* PS PS: Note to myself - make sure everything is reproducible!!!!
* In progress: Testing stability of convergence criteria. Next: Test cases 1) Lipid bilayer 2) Sphere with solutes
* Iterative method works --> Optimizing and validating/testing it remains in progress
* Should use Saga cluster node(s), not Betzy. 
* Hands on work with the code. 
* Note: I will first implement the painted version of q/eps first. When the code works, I can try the q_i/eps_i version --> Discarded this idea. A drawback with painting first, is that it might be inefficient, but it seems "safer" (as of August).
    
** Find format for master thesis (not important atm).
** Make an early disposition for thesis (motivation, existing knowledge , opportunities on what to focus on etc) (not important atm)
** Reading chap 22 in mol driving forces (downprioritize atm, getting hands-on experience first)

Otherwise:
* I will also work a bit on machine learning for a course this fall, so that will take some of my time as well.


## The key steps are:
- [x] 0. Install packages
- [x] 1. Prepare branch and local cloned repository
- [x] 2. Get an overview of the workflow of the algorithms and code
    - [x] understand TomI files
    - [x] load dielectric of different particles to h5
    - [x] read through main.py
    - [x] make an overview over variables in main.py
    - [x]  more detailed description of its main machinery in a flow chart
    - [x] Document functions in field.py
(- [x] Suggest pseudocode for solving the generalized poisson equation (should be close to algo given, so maybe not important))
- [x] 3. Be ready to do implementation by August (optimally)
    - [x] See coding progress down below
- [ ] 4. Optimize and get run results
    - [x] Validate the code
    - [ ] Get error estimates for test (known) case
    - [x]  Run a relevant case and investigate the results

Literature to read
- [ ] Solving generalized poisson equation (not Fourier Space)
- [x] Bore, Cascella (2019) Mesoscale Electrostatics Driving Particle Dynamics in Nonhomogeneous Dielectrics
- [x] Bore, Cascella (2020)Hamiltonian and Alias-Free Hybrid Particle-Field Molecular Dynamics
- [x] Molecular Driving Forces, at least chapter 
    * [x] chap 20; electrostatics
    * [x] chap 21; poisson equation, work and image charges method
- [ ] Bonus chapters (seems interesting)
    * [ ] chap 22; Electrochemical Equilibrum, Nernst equation and Born energy
    * [ ] chap 23, 

## Questions and Answers (Q&A)

| Question   | Answer   | Comment   |
| --------   | --------   | --------    |
|I tried investigating the COULK_GMX = 138.935458 variable, but I haven't found any units where the coloumb constant has this value.| It is in GROMACS units, and each transfer must have a scaling like this, at least wrt. charges| See eg.https://en.wikipedia.org/wiki/Coulomb_constant . In Gaussian units k_e = 1, but also gromacs documentation|
|in main.py: Why integrate twice for velocities in loop for inner rRespa?|Because we do not yet have the velocities due to the bond forces and angle forces, thus we need to integrate velocities from this contrib. before continuing. |text| 
|What is this import line:```from types import ModuleType as moduleobj``` ?| text  |  (Not important) |
| Do we need to change velocity verlet for electrostatics? | Xinmeng:"No, I think now you do not need to worry about it. Tthe electric forces will be added with other forces, and the total force will be used in the Velocity Verlet       | (I know there is a version out there which conserves electrostatic forces)       |
|Do we know that energy is conserved for sure in the existing code? |   Xinmeng:"Yes, Morten has systemtacally tested it. We are sort of confident about this. You will get more details when you use the code" | text | 
|Sigbjørn had some conserns about parallelism and FFT in the iterative method. Input/Thoughts? |Xinmeng "We will discuss about this when you (or we together) write out a psudo code; now you don not need to consider the effciency"| text |




## 0.0 Coding Pogress and thoughts
* PS: The dielectric const is used in field.py through config. This would need to be adressed/changed for the non-constant dielectric --> Here you would solve the poisson equation iteratively instead
* [x] In HyMD-2021 in test-xinmeng-electrostatic in run-dppc: added load-hdf5-add-charge-dielectric-info.py and the dielectric data write to h5.
* [x] Need to add some lines around 341 in main.py, to read off the dielectric, add a dielectric flag
* [x] in main.py add pm.create for needed variables (around line 450)
* [x] in main.py around line 500 add args_in and args_recv for dielectric variable. 
* [x] in field.py - see update_field_force_q for clues - (or main.py?) Need to paint q/eps to grid, and paint eps to grid --> nabla eps/eps. Made a new function update_field_force_q_GPE.
* [x] config.toml changes: see lines 56-58. Need to do simple changes here. Bc. PIC_spectral refers to simple poisson(?) 
* [x] Make a function for the iterative method, place it somewhere convenient - e.g. the newly made function update_field_force_q_GPE.
* [x] Testing the GPE method. It WORKS :smiley: 
    * [x] Fixed overflow with correct scaling in transfer.
    * [x] Check if it is consistent with constant dielectric for a larger mesh e.q. 12,12,12
    * [x] Do I need also to use the COULK_GMX constant in the transfer function to take the gradient? Check this!! 
    * [x] Getting correct output
    * [x] optimizing for mesh size and sigma --> investigate effects of painting eps to grid (Obs: there ARE effects!!).
    * [x] Need to decide upon a convergence criteria --> which has little effects on constant case
* [x] Get main.py with MD to work with dielectric.
    * [x] store-static-with-charge-dielectric
    * [x] Added necessary steps/if tests in outer rrespa steps
* [x] Changed how dielectric is spread to field (grid) paint directly --> weighted mean
* [x] Consider storing the dielectric differently in h5 file for effieciency, if possible.
    * [x] Change the input_parser.py
    * [x] Added Dielectricic_type as dataclass
     in force.py, PSSS: might not need this!!!
     * [x] Dielectric can now be given as a toml dictionary, using dielectric_type = (input dict)
     * [x] check dielectric for error handling?
    * [ ] Can be added --> save to output + print to terminal 
* [ ] Optimizing iterative method.
    * [ ] finding the right factor (1 or 1/c)
    * [ ] Make choice of convergence criteria stable and flexible
* [x] Consider tweaking params for the dielectric. You might need (two different mesh sizes in the toml file, and also) two different sigmas. NOTE: We work with the same sigmas as of now.

* Dielectric paint final decision: 
    *  [x] Made the dielectric into toml convention. 
    * [x] cleanup of different version to paint dielectric.
    * [x] checked that the type is matched with phi
    *  Changes made in  main.py, field.py, input_parser, file_io.py, config.toml and load-h5....py .
    
** [ ] Make separate arrrays for non-zero and zero charges. NOTE: This will be a backup case if needed.

** ([ ] in main.py Around line 600, do layout phi_dielectric ? pm.decompose) I don't think this is needed. At least not atm.

## 0.1 prepare the local DE #Done

> able to run HyMD in local MAC 

```bash 
# Download and install anaconda
  conda create -n py38(88) python=3.8 
 
  Important to do this first and activate

# To activate this environment, use
#     $ conda activate py38
# To deactivate an active environment, use
#     $ conda deactivate
#

install follow the Betzy, when pip3 fails use conda install

however the h5py is not mpi enabled

then 
conda install ( -c conda-forge) "h5py>=2.9=mpi*"
then 
import h5py not found
then 
conda uninstall ( -c conda-forge) "h5py>=2.9=mpi*"
then install 
conda install -c conda-forge "h5py>=2.9=mpi*"
then works 



in dppc folder
bash run.sh 
```

- installed the code in LOCAL ENVIRONMENT  py3888 after many trials and errors.

## 0.2 Overview of algorithm flow

Overview
* Work in the local branch elecGPE
* write in hymd
* Test in test folder
---
How to compile and run code:

```bash
conda activate py3888 # activate local environment
cd [...] test-xinmeng-electrostatic/run-dppc # change directory to correct path
bash run-nocharge.sh ## run a case 
bash run.sh ## run another case
vi run.txt # get more details about the case
vi run-nocharge.sh ## specifies dir and  number of cores
vi run.sh        ## specifies dir number of cores
## you can edit directly in terminal with vi

```
---
### General structure of the code

1) parse inputfiles 
- .tomi
- .h5
-> Read and distribute (parallelism)
2) Setup of arrays
3) Update field
4) Compute field force 
5) Compute band force (spring forces etc)
6) Loop over n time steps
-> velocity verlet
-> compute energy 
-> output .h5  (Not the same as input, we make a new file)

---
Detailed work flow in main.py
![](https://i.imgur.com/jldcigv.png)
![](https://i.imgur.com/X84DvXN.png)

---
in field.py: The idea behind PE solution with *constant* dielectric
![](https://i.imgur.com/XzVGQfr.png)

---
Gromacs unit factor 
![](https://i.imgur.com/0COGYlO.png)

https://en.wikipedia.org/wiki/Gaussian_units

https://en.wikipedia.org/wiki/Relative_permittivity#Terminology

https://manual.gromacs.org/documentation/2019/reference-manual/definitions.html

---
### Important functions, variables and dimensions
(Maybe not that important, can possibly fill out as I go. Also see the work flow of main.py). At least describe what I added myself.

| Variable name | Dimension | Usage/relevant equation |
| -------- | -------- | -------- |
| dielectric     | (20336,) (#particles)     | setting the dielectric (epsilon) for the lipid bilayer and water molecules. Will be used later spread as a density to mesh |
|phi_q_eps|3d mesh - decomp|Non-polarization part of GPE: charge density divided by the dielectric|
|phi_eps|3d mesh - decomp|dielectric (epsilon) spread to mesh as density. For later use in phi_eta|
|phi_eta|3d mesh, for each d a vector of len 3 - decomp|Density of dielectric fraction nabla eps/eps, used in the iterative method and scaling factor of the polarization charge density|
|phi_pol,phi_pol_prev|3d mesh - decomp|Density of the polarization charge used in the iterative method
|update_field force_q_GPE|txt|painting the variables to grid and solving GPE iteratively|



---
A lot of the relevant code for electrostatics can be found in field.py in hymd. Might need to add input on these functions.

---
| Code lines | Purpose/functionality | Functions and contructs used  |
| -------- | -------- | -------- |
| 85     | spread the densities to grid     | pm.paint     |
| 54-89     | get the electrostatic forces  | update_field_force_q    |
| 90-116     | calculate transfer function (see paper)     | phi_transfer_function     |
| 120-196     | Get energies from the field force    | update_field_force_energy_q    |
| update_field     | get phi and update fields    | Text     |
| Text     | Text     | Text     |
| Text     | Text     | Text     |
| Text     | Text     | Text     |


---
### Structure of the toml config files
* structured into meta, particles and simulation...
* need to change (see lines 56-58)
  1) Coulombtype
  2) Dielectric constant

## 0.3 Outline of iterative method
![](https://i.imgur.com/K1Dc8pS.png)
![](https://i.imgur.com/jjpyX5Z.png)
![](https://i.imgur.com/8UCp4Au.png)


## 0.4 Some Literature to electrostatics and Fourier Transform
- :star: https://aip.scitation.org/doi/pdf/10.1063/1.4939125 , Solving generalized poisson equation (not Fourier Space)
    - This article is also relevant with regards to comparing results for the iterative method
- https://pubs.acs.org/doi/10.1021/acs.jctc.8b01201 , Mesoscale Electrostatics Driving Particle Dynamics in Nonhomogeneous Dielectrics
    - I can possibly compare with results from this article. Here they use finite differences to find a solution for a variable dielectric.
- https://aip.scitation.org/doi/full/10.1063/5.0020733 , Hamiltonian and alias-free hybrid particle–field molecular dynamics
- https://physics.stackexchange.com/questions/14639/how-is-the-saddle-point-approximation-used-in-physics, A note about the saddle point approx used for mean field potential $v_{ext}$.

## Links and packages
- Wikipedia has all the answers  :) https://en.wikipedia.org/wiki/Fourier_transform
- the coloumb constant: https://en.wikipedia.org/wiki/Coulomb_constant
- pmesh pm documentation: http://rainwoodman.github.io/pmesh/pmesh.pm.html
- Apply documentation. Seems possible to use a cdot function for inner products, instead of writing my own. Otherwise, see existing code as well. http://rainwoodman.github.io/pmesh/pmesh.pm.html?highlight=apply#pmesh.pm.BaseComplexField.apply
- A page about GROMACS units: https://manual.gromacs.org/documentation/2019/reference-manual/definitions.html 

## 0.6 Macroscopic properties and other observables to test for when the code is done
* Write something here 
* Sigbjørn might have some interesting ideas
* Get close to experiments? 
* Compare with results from other simulation models
* PS: Be critical about other articles and their results as well. 
* Possible shortcomings
* Possible gains or future usage

## Tips
- visualize field 
```python 
import pyvista as pv
data = pv.wrap(_e_potential.value)
outfile = 'xxxx.vtk'
data.save(out_file)

### vtk file can be visualized by Paraview
```
