# 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

# 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:

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

---
in field.py: The idea behind PE solution with *constant* dielectric

---
Gromacs unit factor

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

## 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

import pyvista as pv
data = pv.wrap(_e_potential.value)
outfile = 'xxxx.vtk'
data.save(out_file)

### vtk file can be visualized by Paraview